From statistical modelling to scientific report writing
DSP Handbook for PSB2024
Author
Anna B. Neuheimer
Published
10-11-2024
In this chapter you will learn:
how to report your hypothesis testing results (including visualizations) to communicate what your model says about your hypothesis.
1 Reporting your hypothesis testing results
To report your model, you want to communicate a number of things:
your best-specified model(s)
how well your model explains your response
your modelled effects
In a few weeks, we will discuss communicating your hypothesis testing (including results) more formally - i.e. the text and visuals you will use to write your own scientific reports and papers. In the current chapter, we will gather many of the pieces needed for that task.
2 Introducing the examples:
In this chapter, we will report model results for four example best-specified models that showcase a range of different hypotheses.
If you are re-reading this chapter with your own hypothesis you are testing, look for the example that shares the most structure with your own hypothesis.
The examples are:
Example 1: Resp1 ~ Cat1 + 1
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Example 3: Resp3 ~ Cont4 + 1
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
where
Resp# is your response variable.
Cat# indicates a categorical predictor
Cont# indicates a continuous predictor
Note that each example is testing a different hypothesis by fitting a different model to a different data-set. The four examples are not related to one another.
Note also that the examples will all involve a normal error distribution assumption. I will point out when your process would be different for a different error distribution assumption (e.g. Gamma, poisson, and binomial)
Finally, note that I am not going through all the steps that got us to these example best-specified models, but assume we went through the previous steps (Response, Predictors, Hypothesis, Starting model, Model validation and Hypothesis testing) as we practiced in class to choose each best-specified model.
Example 1: Resp1 ~ Cat1 + 1
Example 1: Resp1 ~ Cat1 + 1
For example #1, assume your best-specified model says that your response variable (Resp1) is explained by a categorical predictor (Cat1):
Resp1 ~ Cat1
Your best-specified model for example #1 is in an object called bestMod1:
bestMod1
Call: glm(formula = Resp1 ~ Cat1 + 1, family = gaussian(link = "identity"),
data = myDat1)
Coefficients:
(Intercept) Cat1Sp2 Cat1Sp3
20.955 5.877 -4.510
Degrees of Freedom: 99 Total (i.e. Null); 97 Residual
Null Deviance: 7083
Residual Deviance: 5210 AIC: 687.1
and the dredge() table you used to pick your bestMod1 is in dredgeOut1:
dredgeOut1
Global model call: glm(formula = Resp1 ~ Cat1 + 1, family = gaussian(link = "identity"),
data = myDat1)
---
Model selection table
(Intrc) Cat1 df logLik AICc delta weight
2 20.95 + 4 -339.547 687.5 0.00 1
1 21.79 2 -354.905 713.9 26.42 0
Models ranked by AICc(x)
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
For Example #2, assume your best-specified model says that your response variable (Resp2) is explained by two categorical predictors (Cat2 & Cat3) as well as the interaction between the predictors:
Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Your best-specified model for example #2 is in an object called bestMod2:
For Example #3, assume your best-specified model says that your response variable (Resp3) is explained by one continuous predictor (Cont4):
Resp3 ~ Cont4 + 1
Your best-specified model for example #3 is in an object called bestMod4:
bestMod3
Call: glm(formula = Resp3 ~ Cont4 + 1, family = gaussian(link = "identity"),
data = myDat3)
Coefficients:
(Intercept) Cont4
-226.9 260.1
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 276100000
Residual Deviance: 64930000 AIC: 1628
that was fit to data in myDat4:
str(myDat3)
'data.frame': 100 obs. of 2 variables:
$ Resp3: num 3668 2774 4594 1094 2100 ...
$ Cont4: num 13.93 12.15 15.4 4.87 6.64 ...
and the dredge() table you used to pick your bestMod3 is in dredgeOut3:
dredgeOut3
Global model call: glm(formula = Resp3 ~ Cont4 + 1, family = gaussian(link = "identity"),
data = myDat3)
---
Model selection table
(Intrc) Cont4 df logLik AICc delta weight
2 -226.9 260.1 3 -811.080 1628.4 0.00 1
1 2358.0 2 -883.457 1771.0 142.63 0
Models ranked by AICc(x)
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
For Example #4, assume your best-specified model says that your response variable (Resp4) is explained by one categorical predictor (Cat4) and one continuous predictor (Cont6) as well as the interaction between the predictors:
Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
Your best-specified model for example #3 is in an object called bestMod4:
This section is similar regardless of your model structure (e.g. error distribution assumption). The examples below all assume a normal error distribution assumption but you can use the same process for models of any structure.
Reporting your best-specified model means reporting the terms - the predictors and any interactions - that are in your model. It is good practice to present the model along with the results from the model selection. In this way, you can include multiple best-specified models if there is evidence that more than one might be useful. Depending on your hypothesis and results, you may want to present all models within ∆AICc < 2 of the best-specified model, or all models with any Akaike weight, or simply all models.
Remember from the Hypothesis Testing chapter that you can also use the output from dredge() function to report evidence for how you picked the best-specified model. When we discuss communicating your hypothesis testing in a few weeks, I’ll include some functions that can help make quick work of making a model selection table.
Let’s take a look at how you do this with our examples:
Example 1: Resp1 ~ Cat1 + 1
Example 1: Resp1 ~ Cat1 + 1
Global model call: glm(formula = Resp1 ~ Cat1 + 1, family = gaussian(link = "identity"),
data = myDat1)
---
Model selection table
(Intrc) Cat1 df logLik AICc delta weight
2 20.95 + 4 -339.547 687.5 0.00 1
1 21.79 2 -354.905 713.9 26.42 0
Models ranked by AICc(x)
For Example 1, you will report that your best-specified model is Resp1 ~ Cat1 + 1, i.e. that there is evidence that Cat1 explains variability in Resp1. This was chosen as the best-specified model as it had the lowest AICc and the next highest rank model had a ∆AICc of 26.42 (i.e. ∆AICc > 2). The best-specified model had an Akaike weight of 1.
For Example 2, you will report that your best-specified model is Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1, i.e. that there is evidence that Cat2 and Cat3 explain variability in Resp2, and that there is an interaction between the effect of Cat2 and Cat3 on your response - i.e. the effect of Cat2 on Resp2 depends on the value of Cat3.
This was chosen as the best-specified model as it had the lowest AICc and the next highest rank model had a ∆AICc of 6.55 (i.e. ∆AICc > 2). The best-specified model had an Akaike weight of 0.964.
Example 3: Resp3 ~ Cont4 + 1
Example 3: Resp3 ~ Cont4 + 1
Global model call: glm(formula = Resp3 ~ Cont4 + 1, family = gaussian(link = "identity"),
data = myDat3)
---
Model selection table
(Intrc) Cont4 df logLik AICc delta weight
2 -226.9 260.1 3 -811.080 1628.4 0.00 1
1 2358.0 2 -883.457 1771.0 142.63 0
Models ranked by AICc(x)
For Example 3, you will report that your best-specified model is Resp3 ~ Cont4 + 1, i.e. that there is evidence that Cont4 explains variability in Resp3.
This was chosen as the best-specified model as it had the lowest AICc and the next highest rank model had a ∆AICc of 142.63 (i.e. ∆AICc > 2). The best-specified model had an Akaike weight of 1.
For Example 4, you will report that your best-specified model is Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1. This model says that there is evidence that Cat5 and Cont6 explain variability in Resp4, and that there is an interaction between Cat5 and Cont6 - i.e. the effect of Cont6 on Resp4 depends on the value of Cat5.
This was chosen as the best-specified model as it had the lowest AICc and the next highest rank model had a ∆AICc of 19.39 (i.e. ∆AICc > 2). The best-specified model had an Akaike weight of 1.
4 Reporting how well your model explains your response
This section is similar regardless of your model structure (e.g. error distribution assumption). The examples below all assume a normal error distribution assumption but you can use the same process for models of any structure.
In this section you will determine
how much deviance in your response is explained by your model
how important each of your predictors is in explaining that deviance
4.1 how much variation in your response is explained by your model
If you will recall, your whole motivation for pursuing statistical modelling was to explain variation in your response. Thus, it is important that you quantify how much variation in your response you are able to explain by your model.
Note that we will discuss this as “explained deviance” rather than “explained variation”. This is because the term “variance” is associated with models where the error distribution assumption is normal, whereas deviance is a more universal term.
When you have a normal error distribution assumption and linear shape assumption, you can capture the amount of explained deviance simply by comparing the variance in your response (i.e. the starting variation before you fit your model - the null variation) vs. the variance in your model residuals (i.e. the remaining variation after you fit your model - the residual variation) as the \(R^2\):
From this equation, you can see how, if your model is able to explain all the variation in the response, the residual variation will be zero, and \(R^2 = 1\). Alternatively, if the model explains no variation in the response the residual variation equals the null variation and \(R^2 = 0\).
For models with other error distribution and shape assumptions, you need another way of estimating the goodness of fit of your model. You can do this through a pseudo \(R^2\).
One useful pseudo \(R^2\) is called the Likelihood Ratio \(R^2\) or \(R^2_{LR}\). The \(R^2_{LR}\) compares the likelihood of your best-specified model to the likelihood of the null model:
where \(n\) is the number of observations, \(log𝓛(model)\) is the log-likelihood of your model, and \(log𝓛(null)\) is the log-likelihood of the null model. The \(R^2_{LR}\) is the type of pseudo \(R^2\) that shows up in your dredge() output when you add extra = "R^2" to the dredge() call. You can calculate \(R^2_{LR}\) by hand, read it from our dredge() output, or estimate it using r.squaredLR() from the MuMIn package:
Note here that two values of \(R^2_{LR}\) are reported. The adjusted pseudo \(R^2\) given here under attr(,"adj.r.squared") has been scaled so that \(R^2_{LR}\) can reach a maximum of 1, similar to a regular \(R^2\)1.
Note that the two estimates are similar: One nice feature of the \(R^2_{LR}\) is that it is equivalent to the regular \(R^2\) when our model assumes a normal error distribution assumption and linear shape assumption, so you can use \(R^2_{LR}\) for any of the models that we are discussing in class.
4.2 how important each of your predictors is in explaining that variation
When you have more than one predictor in your model. you may also want to report how relatively important each predictor is to explaining deviance in your response. This is also called “partitioning the explained deviance” to the predictors.
Aside: what we won’t be doing
To get an estimate of how much deviance in your response one particular predictor explains, you may be tempted to compute the explained deviance (\(R^2\)) estimates of models fit to data with and without that particular predictor. Let’s try with our Example 4:
If you want to get an estimate as to how much response deviance the Cont6 predictor explains, you might compare the \(R^2\) of a model with and without the Cont6 predictor.
Let’s compare comparing model #4 (that includes Cont6) and model #2 (that doesn’t include Cont6):
R2.mod4 <- (dredgeOut4$`R^2`[2]) # model #4 is the second row in dredgeOut4R2.mod2 <- (dredgeOut4$`R^2`[3]) # model #2 is the third row in dredgeOut4diffR2 <- R2.mod4 - R2.mod2 # find estimated contribution of Cont to explained deviancediffR2
[1] 0.02429225
So it looks like 2.4% of the variability in Resp4 is explained by Cont6.
But what if instead you chose to compare model #3 (that includes Cont6) and model #1 (that doesn’t include Cont6):
R2.mod3 <- (dredgeOut4$`R^2`[4]) # model #3 is the fourth row in dredgeOut4R2.mod1 <- (dredgeOut4$`R^2`[5]) # model #1 is the fifth row in dredgeOut4diffR2 <- R2.mod3 - R2.mod1 # find estimated contribution of Cont to explained deviancediffR2
[1] 0.08480009
Now it looks like 8.5% of the variability in Resp4 is explained by Cont6! Quite a different answer! Your estimates of the contribution of Cont6 to explaining the response deviation don’t agree because of collinearity among our predictors2.
There are a few options you can use to get a robust estimate of how much each predictor is contributing to explained deviance in your response.
One option for partitioning the explained deviance when you have collinearity among your predictors is hierarchical partitioning. Hierarchical partitioning estimates the average independent contribution of each predictor to the total explained variance by considering all models in the candidate model set3. This method is beyond the scope of the course but see the rdacca.hp package for an example of how to do this.
Another method (that we will be using) for estimating the importance of each term (predictor or interaction) in your model is by again looking at the candidate model set ranking made by dredge(). Here you can measure the importance of a predictor by summing up the Akaike weights for any model that includes a particular predictor. The Akaike weight of a model compares the likelihood of the model scaled to the total likelihood of all the models in the candidate model set. The sum of Akaike weights for models including a particular predictor tells you how important the predictor is in explaining the deviance in your response. You can calculate the sum of Akaike weights directly with the sw() function in the MuMIn package:
sw(dredgeOut4)
Cat5 Cont6 Cat5:Cont6
Sum of weights: 1 1 1
N containing models: 3 3 1
Here we can see that all model terms (the predictors Cat5 and Cont6 as well as the interaction) are equally important in explaining the deviance in Resp4 (they appear in all models that have any Akaike weight).
Let’s look at these two steps
how much deviance in your response is explained by your model
how important each of your predictors is in explaining that deviance
with our examples:
Example 1: Resp1 ~ Cat1 + 1
Example 1: Resp1 ~ Cat1 + 1
how much deviance in your response is explained by your model
With Example 1, you have only one predictor (Cat1) and so this predictor is responsible for explaining all of the variability in your response (Resp1). You can see that it appears in all models with any weight with your sw() function from the MuMIn package.
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
how much deviance in your response is explained by your model
Cat2 Cat3 Cat2:Cat3
Sum of weights: 1.00 1.00 0.96
N containing models: 3 3 1
Here you can see that Cat2 and Cat3 are equally important in explaining the deviance in Resp2 (they appear in all models that have any Akaike weight), while the interaction term between Cat2 and Cat3 is less important (only appearing in one model with Akaike weight, though this is the top model).
Example 3: Resp3 ~ Cont4 + 1
Example 3: Resp3 ~ Cont4 + 1
how much deviance in your response is explained by your model
R2 <-r.squaredLR(bestMod3)R2
[1] 0.764852
attr(,"adj.r.squared")
[1] 0.764852
The best-specified model explains 76.5% of the deviance in Resp3.
how important each of your predictors is in explaining that deviance
dredgeOut3
Global model call: glm(formula = Resp3 ~ Cont4 + 1, family = gaussian(link = "identity"),
data = myDat3)
---
Model selection table
(Intrc) Cont4 df logLik AICc delta weight
2 -226.9 260.1 3 -811.080 1628.4 0.00 1
1 2358.0 2 -883.457 1771.0 142.63 0
Models ranked by AICc(x)
sw(dredgeOut3)
Cont4
Sum of weights: 1
N containing models: 1
With Example 3, you have only one predictor (Cont6) and so this predictor is responsible for explaining all of the variability in your response (Resp3). You can see that it appears in all models with any weight with your sw() function from the MuMIn package.
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
how much deviance in your response is explained by your model
R2 <-r.squaredLR(bestMod4)R2
[1] 0.855657
attr(,"adj.r.squared")
[1] 0.8567834
The best-specified model explains 85.6% of the deviance in Resp2.
how important each of your predictors is in explaining that deviance
Cat5 Cont6 Cat5:Cont6
Sum of weights: 1 1 1
N containing models: 3 3 1
Here you can see that Cat5, Cont6 and the interaction Cat5:Cont6 are all equally important in explaining the deviance in Resp4 (they appear in all models that have any Akaike weight).
5 Reporting your modelled effects
In the last chapter of these notes, we discussed testing your hypothesis.
Remember that testing your hypothesis is looking for evidence of effects of your predictors on your response. By testing your hypothesis, you determined that there was evidence that the effects of your predictors on your response are non-zero.
Here you will gather the information to communicate visually and quantitatively the magnitude and direction of your effects and determine (in some cases) where effects are similar or different from one another, answering questions like:
how much does a change (effect) in your predictor change your response?
is that change (effect) positive (your response) or negative (your response decreases)?
is that change (effect) similar across levels of your categorical predictor?
You will report your modelled effects in two ways:
Visualizing your model effects: Here, you will make plots of your model effects including the effects themselves, uncertainty around the effects and your observations. These plots will help communicate your modelled effects as well as how well your model fits your observations.
Quantifying your model effects: Here, you will report (in numbers) the magnitude and direction of your modelled effects along with the uncertainty. Exactly how you do that will depend on the structure of your model (e.g. the error distribution assumption). We will go through a number of examples here. You will also look for evidence of whether modelled effects of a categorical predictor differ across all levels.
5.1 Visualizing your model effects
This section is similar regardless of your model structure (e.g. error distribution assumption). The examples below all assume a normal error distribution assumption but you can use the same process for models of any structure.
Here, you will visualize how each predictor is affecting the response by drawing how your modelled fitted values (the estimates of your response made by the model) change with your predictors, along with a measure of uncertainty around those fitted values. You also want to include your data (actual observations of your response and predictors) so you can start to communicate how well your model captures your hypothesis, and what amount of deviance in your response is explained by your model.
There are a lot of different R packages that make it easy to quickly visualize your model. We will focus on two methods that will allow you to make quick visualizations of your model and/or customize the figures as you would like.
visualizing using the visreg package: This will let you get a quick look at your model object with a limited amount of customization.
vizualizing “by hand”: Here, “by hand” is a bit of a silly description as R will be doing the work for you. What I mean by “by hand” is that you will be building the plot yourself by querying your model object. This method is very flexible and will lead you to a fully customizable visualization. This process involves i) choosing the values of your predictors at which to make estimates of your fitted values, ii) using predict() to use your model to estimate your response variable at those values of your predictors, and iii) use the model estimates to plot your model fit. Aside: this is similar to what we will be doing when we learn about model predictions in the next chapter of these notes.
Below are visualizations for each of our examples.
Example 1: Resp1 ~ Cat1 + 1
Example 1: Resp1 ~ Cat1 + 1
Visualizing with the visreg package
library (visreg) # load visreg package for visualizationlibrary(ggplot2) # load ggplot2 package for visualizationvisreg(bestMod1, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cat1", # predictor on x-axis#by = ..., # if you want to include a 2nd predictor plotted as colour#breaks = ..., # if you want to control how the colour predictor is plotted#cond = , # if you want to include a 3rd predictor#overlay = TRUE, # to plot as overlay or panels, when there is #rug = FALSE, # to turn off the rug. The rug shows you where you have positive (appearing on the top axis) and negative (appearing on the bottom axis) residualsgg =TRUE)+# to plot as a ggplot, with ggplot, you will have more control over the look of the plot.ylab("Response, (units)")+# change y-axis labelxlab("Cat1, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
Unfortunately due to limitations of the visreg package, you can not easily add you observations onto plots where the x-axis is a categorical predictor. But that’s ok, because there are other options…
Visualizing by hand
To plot by hand, you will
first make a data frame containing the value of your predictors at which you want to plot effects on the response:
# Set up your predictors for the visualized fitforCat1<-unique(myDat1$Cat1) # every value of your categorical predictor# create a data frame with your predictorsforVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of predictors
Next, you will use the predict() function4 to get the model estimates of your response variable at those values of your predictors:
# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod1, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame
Finally, you will use these model estimates to make your plot:
library(ggplot2)ggplot() +# start ggplotgeom_point(data = myDat1, # add observations to your plotmapping =aes(x = Cat1, y = Resp1), position=position_jitter(width=0.1)) +# control position of data points so they are easier to see on the plotgeom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),linewidth=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat1, y = Fit), shape =22, # shape of pointsize =3, # size of pointfill ="white", # fill colour for plotcol ='black') +# outline colour for plotylab("Resp1, (units)")+# change y-axis labelxlab("Cat1, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Visualizing with the visreg package
As you have more than one predictor in Example 2, there are a number of ways you can visualize your effects.
Here is an example of a visualization with Cat2 on the x-axis and the effects of Cat3 plotted as separate panels:
library (visreg) # load visreg package for visualizationlibrary(ggplot2) # load ggplot2 package for visualizationvisreg(bestMod2, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cat2", # predictor on x-axisby ="Cat3", # if you want to include a 2nd predictor plotted as colour#breaks = ..., # if you want to control how the colour predictor is plotted#cond = , # if you want to include a 3rd predictoroverlay =FALSE, # to plot as overlay or panels, when there is #rug = FALSE, # to turn off the rug. The rug shows you where you have positive (appearing on the top axis) and negative (appearing on the bottom axis) residualsgg =TRUE)+# to plot as a ggplot, with ggplot, you will have more control over the look of the plot.ylab("Response, (units)")+# change y-axis labelxlab("Cat2, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
# tip: we can't include our observations on this plot due to the limitations of the visreg package when plotting with a categorical predictor on the x-axis
And here is an example of a visualization with Cat2 on the x-axis and the effects of Cat3 overlayed on the plot as different colours for each category (level) of Cat3:
visreg(bestMod2, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cat2", # predictor on x-axisby ="Cat3", # if you want to include a 2nd predictor plotted as colour#breaks = ..., # if you want to control how the colour predictor is plotted#cond = , # if you want to include a 3rd predictoroverlay =TRUE, # to plot as overlay or panels, when there is #rug = FALSE, # to turn off the rug. The rug shows you where you have positive (appearing on the top axis) and negative (appearing on the bottom axis) residualsgg =TRUE)+# to plot as a ggplot, with ggplot, you will have more control over the look of the plot.ylab("Response, (units)")+# change y-axis labelxlab("Cat2, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
# tip: we can't include our observations on this plot due to the limitations of the visreg package when plotting with a categorical predictor on the x-axis
Unfortunately due to limitations of the visreg package, you can not easily add you observations onto plots where the x-axis is a categorical predictor. But that’s ok, because there are other options…
Visualizing by hand
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCat2<-unique(myDat2$Cat2) # every level of your categorical predictorforCat3<-unique(myDat2$Cat3) # every level of your categorical predictor# create a data frame with your predictorsforVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of predictors#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod2, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat2, # add observations to your plotmapping =aes(x = Cat2, y = Resp2, col = Cat3), position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) +# control position of data points so they are easier to see on the plotgeom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plotsize=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat2, y = Fit, fill = Cat3), shape =22, # shape of pointsize=3, # size of pointcol ='black', # outline colour for pointposition=position_dodge(width=0.75)) +# control position of data points so they are easier to see on the plotylab("Resp2, (units)")+# change y-axis labelxlab("Cat2, (units)")+# change x-axis labellabs(fill="Cat3, (units)", col="Cat3, (units)") +# change legend titletheme_bw() # change ggplot theme
Example 3: Resp3 ~ Cont4 + 1
Example 3: Resp3 ~ Cont4 + 1
Visualizing with the visreg package
Notice how you can include the gg = TRUE argument to plot this as a ggplot type figure. This allows you to add your data onto the visualization of your model.
library (visreg) # load visreg package for visualizationlibrary(ggplot2) # load ggplot2 package for visualizationvisreg(bestMod3, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cont4", # predictor on x-axis#by = ..., # if you want to include a 2nd predictor plotted as colour#breaks = ..., # if you want to control how the colour predictor is plotted#cond = , # if you want to include a 3rd predictor#overlay = TRUE, # to plot as overlay or panels, when there is #rug = FALSE, # to turn off the rug. The rug shows you where you have positive (appearing on the top axis) and negative (appearing on the bottom axis) residualsgg =TRUE)+# to plot as a ggplot, with ggplot, you will have more control over the look of the plot.geom_point(data = myDat3,mapping =aes(x = Cont4, y = Resp3))+ylab("Response, (units)")+# change y-axis labelxlab("Cont4, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
# Note: you CAN include your observations on your plot with visreg when you have a continuous predictor on the x-axis
Visualizing by hand
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCont4<-seq(from =min(myDat3$Cont4), to =max(myDat3$Cont4), by =1)# a sequence of numbers in your continuous predictor range# create a data frame with your predictorsforVis<-expand.grid(Cont4 = forCont4) # expand.grid() function makes sure you have all combinations of predictors. #### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod3, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat3, # add observations to your plotmapping =aes(x = Cont4, y = Resp3)) +# control position of data points so they are easier to see on the plotgeom_line(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cont4, y = Fit),size =1.2) +# control thickness of linegeom_line(data = forVis, # add uncertainty to your plot (upper line)mapping =aes(x = Cont4, y = Upper),size =0.8, # control thickness of linelinetype =2) +# control style of linegeom_line(data = forVis, # add uncertainty to your plot (lower line)mapping =aes(x = Cont4, y = Lower),size =0.8, # control thickness of linelinetype =2) +# control style of lineylab("Resp3, (units)") +# change y-axis labelxlab("Cont4, (units)") +# change x-axis labeltheme_bw() # change ggplot theme
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
Visualizing with the visreg package
As you will see in @sec_quantifying, there is no evidence of an effect of Cont6 on Resp3 when Cat5 = Farm - i.e. the coefficient (slope) is not different from zero. In such cases, you would typically not include the effect on the plot (red line in the next figure), but I include it here to illustrate the visualization techniques.
A plot with Cont6 on the x-axis:
library (visreg) # load visreg package for visualizationlibrary(ggplot2) # load ggplot2 package for visualizationvisreg(bestMod, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cont6", # predictor on x-axisby ="Cat5", # predictor plotted as colour#breaks = 3, # if you want to control how the colour predictor is plotted#cond = , # if you want to include a 4th predictoroverlay =TRUE, # to plot as overlay or panels rug =FALSE, # to include a ruggg =TRUE)+# to plot as a ggplotgeom_point(data = myDat4, # datamapping =aes(x = Cont6, y = Resp4, col = Cat5))+# add data to your plot#ylim(0, 60)+ # adjust the y-axis unitsylab("Response, (units)")+# change y-axis labelxlab("Cont6, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
A plot with Cat5 on the x-axis: Notice how the continuous predictor is represented by different levels in the colours on the plot. Here you’ve asked for three levels with breaks = 3. Note that you can not include your observations on the visreg plot when the x-axis predictor is a category
visreg(bestMod, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cat5", # predictor on x-axisby ="Cont6", # predictor plotted as colourbreaks =3, # if you want to control how the colour predictor is plotted#cond = , # if you want to include a 4th predictoroverlay =TRUE, # to plot as overlay or panels rug =FALSE, # to include a ruggg =TRUE)+# to plot as a ggplotylab("Response, (units)")+# change y-axis labelxlab("Cat5, (units)")+# change x-axis labellabs(fill="Cont6, (units)", col="Cont6, (units)") +# change legend titletheme_bw() # change ggplot theme
You can also specify at which levels the breaks should occur with the breaks = ... argument. For example, you can ask visreg to plot the modelled effects when Cont6 = 400 and Cont6 = 600:
visreg(bestMod, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cat5", # predictor on x-axisby ="Cont6", # predictor plotted as colourbreaks =c(400,600), # if you want to control how the colour predictor is plotted#cond = , # if you want to include a 4th predictoroverlay =TRUE, # to plot as overlay or panels rug =FALSE, # to include a ruggg =TRUE)+# to plot as a ggplotylab("Response, (units)")+# change y-axis labelxlab("Cat5, (units)")+# change x-axis labellabs(fill="Cont6, (units)", col="Cont6, (units)") +# change legend titletheme_bw() # change ggplot theme
Note that you can not include your observations on the visreg plot when the x-axis predictor is a category.
Visualizing by hand
Here’s an example with Cont6 on the x-axis:
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCat5<-unique(myDat4$Cat5) # all levels of your categorical predictorforCont6<-seq(from =min(myDat4$Cont6), to =max(myDat4$Cont6), by =1)# a sequence of numbers in your continuous predictor range# create a data frame with your predictorsforVis<-expand.grid(Cat5 = forCat5, Cont6 = forCont6) # expand.grid() function makes sure you have all combinations of predictors. #### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod4, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat4, # add observations to your plotmapping =aes(x = Cont6, y = Resp4, col = Cat5)) +# control position of data points so they are easier to see on the plotgeom_line(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cont6, y = Fit, col = Cat5),size =1.2) +# control thickness of linegeom_line(data = forVis, # add uncertainty to your plot (upper line)mapping =aes(x = Cont6, y = Upper, col = Cat5),size =0.4, # control thickness of linelinetype =2) +# control style of linegeom_line(data = forVis, # add uncertainty to your plot (lower line)mapping =aes(x = Cont6, y = Lower, col = Cat5),size =0.4, # control thickness of linelinetype =2) +# control style of lineylab("Resp4, (units)") +# change y-axis labelxlab("Cont6, (units)") +# change x-axis labellabs(fill="Cat5, (units)", col="Cat5, (units)") +# change legend titletheme_bw() # change ggplot theme
Here’s an example with Cat5 on the x-axis:
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCat5<-unique(myDat4$Cat5) # all levels of your categorical predictorforCont6<-c(400, 600) # a sequence of numbers in your continuous predictor range# create a data frame with your predictorsforVis<-expand.grid(Cat5 = forCat5, Cont6 = forCont6) # expand.grid() function makes sure you have all combinations of predictors. #### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod4, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat4, # add observations to your plotmapping =aes(x = Cat5, y = Resp4, col = Cont6), position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) +# control position of data points so they are easier to see on the plotgeom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat5, y = Fit, ymin = Lower, ymax = Upper, col = Cont6),position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plotsize=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat5, y = Fit, fill = Cont6), shape =22, # shape of pointsize=3, # size of pointcol ='black', # outline colour for pointposition=position_dodge(width=0.75)) +# control position of data points so they are easier to see on the plotylab("Resp4, (units)")+# change y-axis labelxlab("Cat5, (units)")+# change x-axis labellabs(fill="Cont6, (units)", col="Cont6, (units)") +# change legend titletheme_bw() # change ggplot theme
Aside: Plotting in 3D
You might notice that many of the plots above are communicating 3-dimensions (one response + two predictors) in a 2-dimensional plot. There are other ways of making 3-dimensional plots in R, e.g. with the visreg package using the visreg2d() function in the visreg package:
visreg2d(fit = bestMod4, # your modelxvar ="Cont6", # what to plot on the x-axisyvar ="Cat5", # what to plot on the y-axisscale ="response") # make sure fitted values (colours) are on the scale of the response
or “by hand” using the geom_raster() function in the ggplot2 package:
# Set up your predictors for the visualized fitforCont6<-seq(from =min(myDat4$Cont6), to =max(myDat4$Cont6), by =0.1) # e.g. a sequence of Cont valuesforCat5<-unique(myDat4$Cat5) # every value of your categorical predictor# create a data frame with your predictorsforVis<-expand.grid(Cont6=forCont6, Cat5=forCat5) # expand.grid() function makes sure you have all combinations of predictors# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod4, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame# create your plotggplot() +# start your ggplotgeom_raster(data = forVis, aes(x = Cont6, y = Cat5, fill = Fit))+# add the 3 dimensions as a rastergeom_point(data = myDat4, aes(x = Cont6, y = Cat5, colour = Resp)) # add your data
EXTRA: As you can see, these plots represent the 3rd dimension by using colour. You can also make actual 3 dimensional plots in R with the plotly package. These plots are interactive which makes them more useful than static 3d plots. Click on the plot and move your mouse to rotate the plot!
library(plotly) # load the plotly packageplot_ly(data = forVis, # the data with your model predictions (made above)x =~Cont6, # what to plot on the x axisy =~Cat5, # what to plot on the y axisz =~Fit, # what to plot on the z axistype="scatter3d", # type of plotmode="markers") %>%# type of plotadd_markers(data = myDat4, x =~Cont6, y =~Cat5, z =~Resp4) # add your data
5.2 Quantifying your model effects
How you quantify your model effects varies with your model structure. If it is your first time reading this, read through the examples that present models assuming a normal error distribution assumption. This will help you understand why we are reporting modelled effects in this way. Then, (if relevant) look at the section on the error distribution assumption you are interested in for your model.
In this section, you will learn how to report your modelled effects in numbers. In general, this section will involve answering:
What are your modelled effects (with uncertainty)?
If you also have a categorical predictor with more than two levels, you will also want to answer:
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Let’s examine each of these in turn:
5.2.1 What are your modelled effects (with uncertainty)?
Recall from the text above that your model effects describe the change in your response given a change in your predictor.
When you have a continuous predictor, this effect is the the amount of change in the response that is caused by a unit change in your continuous predictor. Information about this effect is captured in the coefficient (slope) of the linear fit. Expand the text box below to see why.
the math
Here’s a model equation for a generalized linear model:
where * \(g()\) is the link function * \(\mu_i\) is the mean fitted value \(\mu_i\) dependent on \(Pred_i\) * \(\beta_1\) is the coefficient of \(Pred_i\) * \(Pred_i\) is the value of your predictor * \(\beta_0\) is the intercept * \(Resp_i\) is the value of the response variable * \(F\) represents some error distribution assumption
As above, the effect of a continuous predictor on your response is the change in your fitted value (\(g(\mu_i)\)) for a unit change in your predictor (\(Pred_i\)):
\[
\begin{align}
g(\mu|Pred+1) - g(\mu|Pred) &= (\beta_1 \cdot (Pred+1) + \beta_0) - (\beta_1 \cdot Pred_i + \beta_0) \\
&=\beta_1 \cdot Pred + \beta_1 + \beta_0 - \beta_1 \cdot Pred_i - \beta_0 \\
&=\beta_1
\end{align}
\] For linear models with a normal error distribution assumption (i.e. when you are using link = "identity"), the link function \(g()\) is not there. For this reason, you can read the effects of your model directly from the model coefficients (\(\beta_1\), \(\beta_0\)) when you have a normal error distribution assumption. And for this reason, you will need to convert your model coefficients to quantify your modelled effects when you have an error distribution assumption that is not normal/uses a link function that isn’t “identity”.
This is because of the difference between the “link” scale and “response” scale for generalized linear models.
The link scale is where the model fitting happens, and where the evidence for your modelled effects is gathered.
The response scale is back in the real world - these are the actual effects in your response you would expect to see with a change in your predictor.
Recall that you can see the default link functions associated with each error distribution function with ?family:
When you have a categorical predictor, the coefficient represents the intercept of the linear fit and how that intercept changes as you move from one category to another. So if you have a categorical predictor with 3 categories (levels), you will have three coefficients that describe how the mean of your response changes as you move from category to category.
When you have an interaction among your predictors, you will have one coefficient for each combination in your interaction.
Finally, you need to report the uncertainty around your modelled effects so your peers know how much evidence there is for these modelled effects. There are a number of ways to do this. Two that we will use are i) reporting standard errors (SE) which are a measure of uncertainty of the average modelled effect, and ii) the 95% confidence intervals around your modelled effects, which tell you the range your modelled effects would be expected to fall into (with a 95% probability) if you redid your experiment.
These principles are true regardless of your model structure BUT our interpretation of what the effects mean in the real world will vary based on e.g. your error distribution assumptions (see the text box above for an explanation). In the examples below, you will learn how to quantify your modelled effects and uncertainty estimates for models with a normal error distribution assumption. There are sections below to show you how to do the same for models with other error distribution assumptions.
For a quick summary, here is how the method depends on your modelled error distribution assumption:
5.2.2 Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors)
When you have a categorical predictor in your best-specified model, you know there is evidence of an effect for at least one of your levels (categories). This step will tell you where you have evidence that the effects among your factor levels differ. The processes outlined below are the same for any model structure.
5.2.3 Models with a normal error distribution assumption:
Example 1: Resp1 ~ Cat1 + 1
Example 1: Resp1 ~ Cat1 + 1
Recall that Example 1 contains only one predictor and it is categorical:
Resp1 ~ Cat1 + 1
What are your modelled effects (with uncertainty)?
In the chapter on your Starting Model, you found your coefficients in the “Estimate” column of the summary() output of your model:
coef(summary(bestMod1)) # extract the coefficients from summary()
Interpreting the coefficients from this output takes practice - especially for categorical predictors because of the way that R treats categorical predictors in regression. With R’s “dummy coding”, one level of the predictor (here Sp1) is incorporated into the intercept estimate (21) as the reference level. The other coefficients in the Estimate column represent the change in modelled response when you move from the reference level (Sp1) to another level.
The modelled response when Cat1 = Sp2 is (5.9) units higher than this reference level = 21 + 5.9 = 11.8 units.
The modelled Resp1 when Cat1 = Sp3 is -4.5 units lower than the reference level = 21 + 5.9 = 1.4 units.5
So all the information we need is in this summary() output, but not easy to see immediately. An easier way is to use the emmeans6 package which helps you by reporting the coefficients directly. In our case, you use the emmeans package to get the mean value of the response at each level of the predictor. For categorical predictors, you do this with the emmeans() function:
library(emmeans) # load the emmeans packageemmOut <-emmeans(object = bestMod1, # your modelspecs =~ Cat1, # your predictorstype ="response") # report coefficients on the response scaleemmOut
So now you have a modelled value of your response for each level of your categorical predictor - this captures the effect of your categorical predictor on your response. You also need to report uncertainty around this effect, and you can do this by reporting the standard error (SE) also reported by the emmeans() function.
When Cat1 is Sp1, Resp1 is estimated to be 21 ± 1.3 units.
When Cat1 is Sp2, Resp1 is estimated to be 26.8 ± 1.2 units.
When Cat1 is Sp3, Resp1 is estimated to be 16.4 ± 1.3 units.
Note that this is the same information in the summary() output just easier to read.7
Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if I observed a new Resp1 at a particular Cat1, there would be a 95% chance it would lie between the bounds of the confidence limits.
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Since you have a categorical predictor with more than two levels, you will ask if all levels of Cat1 affect Resp1 in the same way - i.e. are the coefficients across factor levels significantly different from one another?
To get evidence about how each level affects your response, you need to test which effects differ from one another among the categories (levels) of your categorical predictor. This is done using multiple comparisons, i.e. you will compare the modelled effect of each level of your categorical predictor vs. each other level of your categorical predictor to determine which are different from each other (called pairwise testing).
This is done by testing the null hypothesis that the modelled effects of each pair are similar to one another, typically rejecting the null hypothesis when P < 0.05. Remember that the P-value is the probability that you observe a value at least as big as the one you observed even though our null hypothesis is true. In this case, you are looking at the value you are observing is the difference between coefficients estimated for two levels of your predictor. The P-value is the probability of getting a difference at least as big as the one you observed even though there is actually no difference between the coefficients (the null hypothesis is true).
A couple of things to note about multiple comparison testing:
Multiple comparison testing on a categorical predictor should only be done after your hypothesis testing has given you evidence that the predictor has an effect on your response. That is, you should only do a multiple comparison test on a predictor if that predictor is in your best-specified model. As this is a test done after your hypothesis testing, it is called a post-hoc8 test.
Multiple comparison testing can be a problem because you are essentially repeating a hypothesis test many times on the same data (i.e. are the effects of Sp1 different than those of Sp2? are the effects of Sp1 different than those of Sp3? are the effects of Sp2 different than those of Sp3?…). These repeated tests mean there is a high chance that you will find a difference in one test purely due to random chance, vs. due to there being an actual difference. To account for this, the multiple comparison tests you will perform have been formulated to correct for this increased error.
Multiple comparison testing is very simple with the emmeans package. It just requires you to hand the output from emmeans() to a function called pairs():
forComp <-pairs(emmOut)forComp
contrast estimate SE df t.ratio p.value
Sp1 - Sp2 -5.88 1.77 97 -3.313 0.0037
Sp1 - Sp3 4.51 1.86 97 2.423 0.0451
Sp2 - Sp3 10.39 1.77 97 5.856 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates
The output shows the results of the multiple comparison (pairwise) testing. The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor. For example, for Sp1 vs. Sp2, there is a 0.37% (p = 0.0037) probability of getting a difference in coefficients at least as big as 26.8 - 21 = r 5.8 even though the null hypothesis (no difference) is true. This value 0.37% (P = 0.0037) is small enough that you can say you have evidence that the effect of Cat1 on Resp1 is different for these two different levels (Sp1andSp2`).
Based on a threshold p-value of 0.05, we can see that:
There is evidence that the value of Resp1 when Cat1 is Sp1 is different (lower) than that when Cat1 is Sp2 as P = 0.0037 is less than P = 0.05.
There is a little evidence that the value of Resp1 when Cat1 is Sp1 is different (higher) than that when Cat1 is Sp3 as P = 0.0451 is less than P = 0.05 (but it is close!). There is some evidence that the value of Resp1 when Cat1 is Sp2 is different (higher) than that when Cat1 is Sp3 as P < 0.00019 is smaller than P = 0.05.10.
Note that the p-values have been adjusted via the Tukey method which adjusts the difference that the two coefficients need to have to allow for the fact that we are making multiple comparisons.11
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Note that the difference between Resp1 when Cat1 is Sp1 vs. Sp3 is so close to overlapping zero. This is reflected in the high P value.
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Recall that Example 2 contains two predictors and both are categorical:
Resp2 ~ Cat2 + Cat3 + Cat2:Cat3 + 1
What are your modelled effects (with uncertainty)?
Again, for categorical predictors, there is a coefficient estimated for each level of the predictor. If there is one or more interactions among predictors, there will be one coefficient for each combination of levels among predictors. Let’s look at the summary() output of your model to understand better:
coef(summary(bestMod2)) # extract the coefficients from summary()
Comparing this output to that of Example 1 shows many more estimated coefficients for Example 2. This is because we have one coefficient for each level of each predictor, as well as coefficients for each combination (the interaction term) of levels of the predictors.
Again, it takes a little practice to read the coefficients from the summary() output. For Example 2:
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 units (the intercept).
The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 365 + 37 = 402 units.
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 365 - 44 = 321 units.
The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 365 - 76 - 44 + 74 = 319 units.
As above, we can use the emmeans package to more easily see these coefficients:
emmOut <-emmeans(object = bestMod2, # your modelspecs =~ Cat2 + Cat3 + Cat2:Cat3, # your predictorstype ="response") # report coefficients on the response scaleemmOut
This is a hard table to navigate as every possible combination of the levels across predictors is compared. This may not be what you want. Instead, you can look for differences of one predictor based on a value of another. For example:
forComp <-pairs(emmOut, # emmeans outputsimple ="Cat2") # contrasts within this categorical predictorforComp
Here you can more easily see the contrasts, and the adjustment to the P-value will be limited to what you actually need. Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other.
For example, the coefficient (predicted Resp2) when Cat3 = Treat1 and Cat2 = TypeA is 365 and this is 37 lower than the coefficient when Cat3 = Treat1 and Cat2 = TypeB, but there is no evidence that the difference has any meaning as P = 0.67 for this comparison.
On the other hand, some combinations of your predictors are statistically different from each other. For example, a comparison of modelled Resp2 when Cat3 = Control and Cat2 = TypeB (548) vs. Cat3 = Control and Cat2 = TypeD (215) shows that they differ by 333 and that there is strong evidence that this difference is meaningful as P < 0.001.
Note that you could also organize this with contrasts by Cat3 instead:
forComp <-pairs(emmOut, # emmeans outputsimple ="Cat3") # contrasts within this categorical predictorforComp
Cat2 = TypeA:
contrast estimate SE df t.ratio p.value
Treat1 - Treat2 43.630 36.4 88 1.198 0.4574
Treat1 - Control -64.116 30.3 88 -2.116 0.0924
Treat2 - Control -107.745 31.9 88 -3.374 0.0031
Cat2 = TypeB:
contrast estimate SE df t.ratio p.value
Treat1 - Treat2 0.921 39.3 88 0.023 0.9997
Treat1 - Control -145.366 30.9 88 -4.711 <.0001
Treat2 - Control -146.287 39.3 88 -3.719 0.0010
Cat2 = TypeC:
contrast estimate SE df t.ratio p.value
Treat1 - Treat2 -29.979 40.1 88 -0.748 0.7357
Treat1 - Control -45.038 28.1 88 -1.605 0.2489
Treat2 - Control -15.060 36.4 88 -0.414 0.9099
Cat2 = TypeD:
contrast estimate SE df t.ratio p.value
Treat1 - Treat2 -22.391 41.4 88 -0.541 0.8514
Treat1 - Control 7.609 34.8 88 0.218 0.9741
Treat2 - Control 30.000 34.8 88 0.861 0.6661
P value adjustment: tukey method for comparing a family of 3 estimates
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Example 3: Resp3 ~ Cont4 + 1
Example 3: Resp3 ~ Cont4 + 1
Recall that Example 3 contains one predictor and the predictor is continuous:
Resp3 ~ Cont4 + 1
What are your modelled effects (with uncertainty)?
For a continuous predictor, the effect is captured in one coefficient that describes the change in the response for a unit change in the continuous predictor. For models with a normal error distribution assumption, this is communicated as the slope.
Let’s look at the summary() output for Example 3:
coef(summary(bestMod3)) # extract the coefficients from summary()
Estimate Std. Error t value Pr(>|t|)
(Intercept) -226.9131 166.09360 -1.366176 1.750105e-01
Cont4 260.0574 14.56593 17.853818 1.438227e-32
This tells you that for a unit change in Cont4, you get a 260 ± 1512.
Note that this same information is in the emmeans package. Because Cont4 is a continuous predictor, you need the emtrends() function, which provides the 95% confidence interval on the estimate:
trendsOut <-emtrends(bestMod3, # your modelspecs =~ Cont4 +1, # your predictorsvar ="Cont4") # your continuous predictorstrendsOut
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
This question is only applies to categorical predictors of which there are none in Example 3.
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
Example 4: Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
Recall that Example 4 contains two predictors, one is categorical and one is continuous:
Resp4 ~ Cat5 + Cont6 + Cat5:Cont6 + 1
What are your modelled effects (with uncertainty)?
The effects estimated for this model will include values of the response (Resp4) at each level of the categorical predictor (Cat5) as well as slopes describing the change in the response when the continuous predictor (Cont6) changes, and a different slope will be estimated for each level in Cat5 (this is because of the interaction in the model).
As above, you can take a look at your modelled coefficients with the summary() output:
The model prediction of Resp4 when Cat5 is Farm and Cont6 is 0 is 98.34 units13 (the intercept).
The model prediction of Resp4 when Cat5 is Urban and Cont6 is 0 is 98.34 - 0.24 = 98.1 units.
The model prediction of Resp4 when Cat5 is Wild and Cont6 is 0 is 98.34 - 0.87 = 97.5 units.
The slope of the relationship between Cont6 and Resp4 when Cat5 is Farm is -0.0029.14
The slope of the relationship between Cont6 and Resp4 when Cat5 is Urban is -0.0029 + 0.015 = 0.0121.
The slope of the relationship between Cont6 and Resp4 when Cat5 is Wild is -0.0029 + 0.031 = 0.0281.
Again, interpreting the coefficients from the summary() output is tedious and not necessary: You can use the emmeans package to give you the modelled response for each level of the categorical predictor (Cat5) directly:
emmOut <-emmeans(object = bestMod4, # your modelspecs =~ Cat5 + Cont6 + Cat5:Cont6 +1, # your predictorstype ="response") # report coefficients on the response scaleemmOut
Note that emmeans() sets our continuous predictor (Cont6) to the mean value of Cont6 (510 units). We can also control this with the at = argument:
emmOut <-emmeans(object = bestMod4, # your modelspecs =~ Cat5 + Cont6 + Cat5:Cont6 +1, # your predictorstype ="response", # report coefficients on the response scaleat =list(Cont6 =0)) # control the value of your continuous predictor at which to make the coefficient estimatesemmOut
By setting at = 0, you get the intercept - i.e. the modelled Resp4 when Cont6 = 0 for each level of Cat5, and this is what is reported in the summary() output.
Similarly, you can get the emmeans package to give you the slope coefficients for the continuous predictor (Cont6) using the emtrends() function, rather than interpreting them from the summary() output:
trendsOut <-emtrends(bestMod4, # your modelspecs =~ Cat5 + Cont6 + Cat5:Cont6 +1, # your predictorsvar ="Cont6") # your continuous predictorstrendsOut
there is no evidence on an effect on Resp4 of Cont6 when Cat5 = Farm (p = 0.44).
there is evidence for a positive effect of Cont6 on Resp4 when Cat5 = Urban. The effect is 0.01215 with a 95% confidence interval of 0.0039 - 0.020516.
there is evidence for a positive effect of Cont6 on Resp4 when Cat5 = Wild. The effect is 0.02817 with a 95% confidence interval of 0.0182 - 0.037718.
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Since you have a categorical predictor with more than two levels, you will want to report where the effects among levels differ from one another.
Note that the effects in Example 4 include slope coefficients (associated with the continuous predictor) and intercept coefficients (associated with the categorical predictor). Because the slope and intercept are dependent on one another (when the slope changes, the intercept will naturally change), you need to first check if the slope coefficients vary across your categorical levels. Only if they don’t vary (i.e. the modelled effects produce lines that are parallel across your categorical levels) would you go on to test if the intercept coefficiens vary across your categorical levels.
You can find out which slopes (i.e. the effect of Cont6 on Resp4) are different across the levels of Cat5 using emtrends() with a request for pairwise testing:
forCompSlope <-pairs(trendsOut)forCompSlope
contrast estimate SE df t.ratio p.value
Farm Cont6509.767 - Urban Cont6509.767 -0.0151 0.00561 94 -2.694 0.0226
Farm Cont6509.767 - Wild Cont6509.767 -0.0309 0.00618 94 -4.991 <.0001
Urban Cont6509.767 - Wild Cont6509.767 -0.0157 0.00646 94 -2.437 0.0437
P value adjustment: tukey method for comparing a family of 3 estimates
Here you see that the slope estimates (the effect of Cont6 on Resp4) for all levels of Cat5 are significantly different from one another (0.0001 ≤ P ≤ 0.044).
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Again: only if all the slopes were similar, would you want to test if the levels of your categorical predictor (Cat4) have significantly different coefficients (intercepts) from each other. You would do this with pairs() on the output from emmeans() (the emmOut object).
5.2.4 Models with a poisson error distribution assumption:
Background
For models with a poisson error distribution assumption, you will typically start with a link function that is the natural logarithm or link = "log". You need to take this link function into consideration when you report your modelled effects.
If you are using link = "log" (as with a poisson error distribution, and sometimes also a Gamma error distribution - more on this to come), you get your coefficient on the response scale by taking e to the coefficient on the link scale (with the R function exp()). This coefficient is called the rate ratio19 and it tells you the % change in the response for a unit change in your predictor.
For a continuous predictor, the effect is captured in one coefficient that describes the change in the response for a unit change in the continuous predictor.
For models with a normal error distribution assumption, this was simply a slope.
For models using a log link (including many models with a poisson error distribution assumption), this is communicated as the rate ratio or the % change in the response for a unit change in your predictor. This allows you to communicate the effect of the continuous predictor while accounting for the curve in the relationship (see visualization of model effects above). This curve occurs because the model is not allowed to go below zero to be consistent with a poisson error distribution assumption.
the math
But why do we convert coefficients from the link to response scale with \(e^x\)20 when link = "log"?
A reminder that we can present our linear model like this:
where * \(g()\) is the link function * \(\mu_i\) is the mean fitted value \(\mu_i\) dependent on \(Pred_i\) * \(\beta_1\) is the coefficient of \(Pred_i\) * \(Pred_i\) is the value of your predictor * \(\beta_0\) is the intercept * \(Resp_i\) is the value of the response variable * \(F\) represents some error distribution assumption
If you have an error distribution assumption (\(F\)), the link function \(g()\) is \(log_e\)21:
The emmeans() function will convert the coefficients (intercepts) from the link to the response scale. You can ask for this with the type = "response" argument.
In contrast, the emtrends() function does not convert the coefficients (slopes) to represent effects on the response scale. This is because emtrends() is reporting the slope of a straight line - the trend line on the link scale. But the line isn’t straight on the response scale.
We need to convert from the link to the response by hand when we use emtrends()^[remember, the coefficients for categorical predictors are also being converted from the link to the response scale. It is just that R is doing it for us in the emmeans() function when we add the argument `type = “response”). How we make the conversion, and interpret the result, will depend on our error distribution assumption:
Here we will go through examples of reporting for models with a poisson error distribution assumption. If you want more details on why you are using the methods below, check the section on “Models with a normal error distribution assumption” (Section 5.2.3) for context.
Examples
Here I have made new examples where the error distribution assumption is poisson. This changes the response variables, indicated by Resp1.p, Resp2.p, etc. to be integers ≥ 0.
Example 1: Resp1.p ~ Cat1 + 1
Example 1: Resp1.p ~ Cat1 + 1
Recall that Example 1 contains only one predictor and it is categorical:
Resp1.p ~ Cat1 + 1
The best-specified model (now with poisson error distribution assumption) is:
bestMod1.p
Call: glm(formula = Resp1.p ~ Cat1 + 1, family = poisson(link = "log"),
data = Dat)
Coefficients:
(Intercept) Cat1StB Cat1StC
6.15524 -0.04869 0.09948
Degrees of Freedom: 49 Total (i.e. Null); 47 Residual
Null Deviance: 139
Residual Deviance: 40.01 AIC: 446
and the dredge() table you used to pick your bestMod1.p is in dredgeOut1.p
dredgeOut1.p
Global model call: glm(formula = Resp1.p ~ Cat1 + 1, family = poisson(link = "log"),
data = Dat)
---
Model selection table
(Intrc) Cat1 df logLik AICc delta weight
2 6.155 + 3 -219.983 446.5 0.00 1
1 6.164 1 -269.477 541.0 94.55 0
Models ranked by AICc(x)
And here’s a visualization of the model effects:
# Set up your predictors for the visualized fitforCat1<-unique(myDat1.p$Cat1) # every value of your categorical predictor# create a data frame with your predictorsforVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of predictors# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod1.p, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data framelibrary(ggplot2)ggplot() +# start ggplotgeom_point(data = myDat1.p, # add observations to your plotmapping =aes(x = Cat1, y = Resp1.p), position=position_jitter(width=0.1)) +# control position of data points so they are easier to see on the plotgeom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),linewidth=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat1, y = Fit), shape =22, # shape of pointsize =3, # size of pointfill ="white", # fill colour for plotcol ='black') +# outline colour for plotylab("Resp1.p, (units)")+# change y-axis labelxlab("Cat1, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
For categorical predictors, you can get the coefficients on the response scale directly with the emmeans() function22 when you specify type = “response”:
library(emmeans) # load the emmeans packageemmOut <-emmeans(object = bestMod1.p, # your modelspecs =~ Cat1 +1, # your predictorstype ="response") # report coefficients on the response scaleemmOut
Cat1 rate SE df asymp.LCL asymp.UCL
StA 471 6.54 Inf 459 484
StB 449 4.32 Inf 440 457
StC 520 5.89 Inf 509 532
Confidence level used: 0.95
Intervals are back-transformed from the log scale
Notice that the column reported is “rate”: these are the coefficients converted back onto the response scale.
Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if I observed a new Resp1.p at a particular Cat1, there would be a 95% chance it would lie between the bounds of the confidence limits.
From this output you see:
When Cat1 is SpA, Resp1.p is estimated to be 471 (95% confidence interval: 459 to 484 units).
When Cat1 is SpB, Resp1.p is estimated to be 449 (95% confidence interval: 440 to 457 units). When Cat1 is SpC, Resp1.p is estimated to be 520 (95% confidence interval: 509 to 532 units).
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Here you can compare the effects of the levels of Cat1 on Resp1.p to see if the effects differ between the levels.
forComp <-pairs(emmOut)forComp
contrast ratio SE df null z.ratio p.value
StA / StB 1.050 0.0177 Inf 1 2.880 0.0111
StA / StC 0.905 0.0162 Inf 1 -5.552 <.0001
StB / StC 0.862 0.0128 Inf 1 -9.968 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates
Tests are performed on the log scale
The output shows the results of the multiple comparison (pairwise) testing.
The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor.
For example:
there is evidence that Resp1.p is significantly larger when Cat 1 = SpA than when Cat1 = SpB (p = 0.0111). Resp1.p is expected to be ~ 5 ± 2% bigger when Cat 1 = SpA than when Cat 1 = SpB.
there is evidence that Resp1.p is significantly smaller when Cat 1 = SpA than when Cat1 = SpC (p = 0). Resp1.p is expected to be ~ 9 ± 2% smaller when Cat 1 = SpA than when Cat 1 = SpC.
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Example 2: Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Example 2: Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1
Recall that Example 2 contains two predictors and both are categorical:
Resp2.p ~ Cat2 + Cat3 + Cat2:Cat3 + 1
The best-specified model (now with poisson error distribution assumption) is:
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCat2<-unique(myDat2.p$Cat2) # every level of your categorical predictorforCat3<-unique(myDat2.p$Cat3) # every level of your categorical predictor# create a data frame with your predictorsforVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of predictors#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod2.p, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat2.p, # add observations to your plotmapping =aes(x = Cat2, y = Resp2.p, col = Cat3), position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) +# control position of data points so they are easier to see on the plotgeom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plotsize=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat2, y = Fit, fill = Cat3), shape =22, # shape of pointsize=3, # size of pointcol ='black', # outline colour for pointposition=position_dodge(width=0.75)) +# control position of data points so they are easier to see on the plotylab("Resp2.p, (units)")+# change y-axis labelxlab("Cat2, (units)")+# change x-axis labellabs(fill="Cat3, (units)", col="Cat3, (units)") +# change legend titletheme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
As above, we can use the emmeans package to see the effects of the predictors on the scale of the response:
emmOut <-emmeans(object = bestMod2.p, # your modelspecs =~ Cat2 + Cat3 + Cat2:Cat3, # your predictorstype ="response") # report coefficients on the response scaleemmOut
Cat2 Cat3 rate SE df asymp.LCL asymp.UCL
TypeA Control 15.20 1.740 Inf 12.14 19.03
TypeB Control 10.86 1.250 Inf 8.67 13.59
TypeC Control 9.67 1.800 Inf 6.72 13.91
TypeD Control 4.25 1.030 Inf 2.64 6.84
TypeA Treat 10.38 0.894 Inf 8.77 12.29
TypeB Treat 4.00 0.707 Inf 2.83 5.66
TypeC Treat 3.62 0.673 Inf 2.52 5.22
TypeD Treat 4.50 1.500 Inf 2.34 8.65
Confidence level used: 0.95
Intervals are back-transformed from the log scale
So now we have a modelled value of our response for each level of our categorical predictors. For example:
The modelled prediction for Resp2.p when Cat2 is TypeA and Cat3 is Treat is 15.2 (95% confidence interval: 12.1 - 19.0 units).
The modelled prediction for Resp2.p when Cat2 is TypeC and Cat3 is Control is 9.7 units (95% confidence interval: 6.7 - 13.9 units).
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
As with Example 1, you can also find out which combinations of predictor levels are leading to statistically different model predictions in Resp2.p:
forComp <-pairs(emmOut, # emmeans outputsimple ="Cat2") # contrasts within this categorical predictorforComp
Here you can more easily see the contrasts. Contrast ratios similar to 1 mean that the predictor levels lead to a similar value of the response.
Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other. For example, the coefficient (predicted Resp2.p) when Cat3 = Control and Cat2 = TypeA (15.2) is 1.4 times (or 40%) more than the predicted Resp2.p when Cat3 = Control and Cat2 = TypeB (10.9), but there is no evidence that the difference has any meaning as P = 0.16 for this comparison.
On the other hand, some combinations of your predictors are statistically different from each other. For example, a comparison of modelled Resp2.p when Cat3 = Control and Cat2 = TypeB (10.9) is 2.6 times (or 160%) higher than when Cat3 = Control and Cat2 = TypeD (4.2) and there is strong evidence that this difference is meaningful as P < 0.0027.
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Example 3: Resp3.p ~ Cont4 + 1
Example 3: Resp3.p ~ Cont4 + 1
Recall that Example 3 contains one predictor and the predictor is continuous:
Resp3 ~ Cont4 + 1
The best-specified model (now with poisson error distribution assumption) is:
bestMod3.p
Call: glm(formula = Resp ~ Cont4 + 1, family = poisson(link = "log"),
data = myDat3.p)
Coefficients:
(Intercept) Cont4
1.577204 0.006674
Degrees of Freedom: 49 Total (i.e. Null); 48 Residual
Null Deviance: 74.74
Residual Deviance: 62.96 AIC: 256.7
that was fit to data in myDat3.p:
str(myDat3.p)
'data.frame': 50 obs. of 2 variables:
$ Cont4 : num 67.5 61.9 71.5 116.6 55.9 ...
$ Resp3.p: int 12 8 4 12 0 2 4 10 8 7 ...
and the dredge() table you used to pick your bestMod3.p is
dredgeOut3.p
Global model call: glm(formula = Resp ~ Cont4 + 1, family = poisson(link = "log"),
data = myDat3.p)
---
Model selection table
(Intrc) Cont4 df logLik AICc delta weight
2 1.577 0.006674 2 -126.339 256.9 0.00 0.992
1 2.069 1 -132.230 266.5 9.61 0.008
Models ranked by AICc(x)
And here’s a visualization of the model effects:
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCont4<-seq(from =min(myDat3.p$Cont4), to =max(myDat3.p$Cont4), by =1)# a sequence of numbers in your continuous predictor range# create a data frame with your predictorsforVis<-expand.grid(Cont4 = forCont4) # expand.grid() function makes sure you have all combinations of predictors. #### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod3.p, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat3.p, # add observations to your plotmapping =aes(x = Cont4, y = Resp3.p)) +# control position of data points so they are easier to see on the plotgeom_line(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cont4, y = Fit),size =1.2) +# control thickness of linegeom_line(data = forVis, # add uncertainty to your plot (upper line)mapping =aes(x = Cont4, y = Upper),size =0.8, # control thickness of linelinetype =2) +# control style of linegeom_line(data = forVis, # add uncertainty to your plot (lower line)mapping =aes(x = Cont4, y = Lower),size =0.8, # control thickness of linelinetype =2) +# control style of lineylab("Resp3.p, (units)") +# change y-axis labelxlab("Cont4, (units)") +# change x-axis labeltheme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
trendsOut <-emtrends(bestMod3.p,specs =~ Cont4 +1, # your predictorsvar ="Cont4") # your predictor of interesttrendsOut
Note that these are the same estimates as we find in the summary() of your model object
coef(summary(bestMod3.p))
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.577203770 0.155306248 10.155443 3.133832e-24
Cont4 0.006674037 0.001934204 3.450534 5.594789e-04
This is an estimate of the modelled effects in the linked world. We need to consider our log link function to report the modelled effects in the response world.
Since you have a poisson error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with:
trendsOut.df <-data.frame(trendsOut) # convert trendsOut to a dataframeRR <-exp(trendsOut.df$Cont4.trend) # get the coefficient on the response scaleRR
[1] 1.006696
We can also use the same conversion on the confidence limits of the modelled effect, for the upper:
RR.up <-exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scaleRR.up
[1] 1.01052
and lower 95% confidence level:
RR.low <-exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scaleRR.low
[1] 1.002887
This tells you that for a unit change in Cont4, Resp3 increases by 1.007x (95% confidence interval is 1.003 to 1.011) which is an increase of 0.67%.
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
This question is only applies to categorical predictors of which their on none in Example 3.
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCat5<-unique(myDat4.p$Cat5) # all levels of your categorical predictorforCont6<-seq(from =min(myDat4.p$Cont6), to =max(myDat4.p$Cont6), by =1)# a sequence of numbers in your continuous predictor range# create a data frame with your predictorsforVis<-expand.grid(Cat5 = Cat5, Cont6 = forCont6) # expand.grid() function makes sure you have all combinations of predictors. #### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod4.p, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat4.p, # add observations to your plotmapping =aes(x = Cont6, y = Resp4.p, col = Cat5)) +# control position of data points so they are easier to see on the plotgeom_line(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cont6, y = Fit, col = Cat5),size =1.2) +# control thickness of linegeom_line(data = forVis, # add uncertainty to your plot (upper line)mapping =aes(x = Cont6, y = Upper, col = Cat5),size =0.4, # control thickness of linelinetype =2) +# control style of linegeom_line(data = forVis, # add uncertainty to your plot (lower line)mapping =aes(x = Cont6, y = Lower, col = Cat5),size =0.4, # control thickness of linelinetype =2) +# control style of lineylab("Resp4, (units)") +# change y-axis labelxlab("Cont6, (units)") +# change x-axis labellabs(fill="Cat5, (units)", col="Cat5, (units)") +# change legend titletheme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
The process for estimating effects of categorical vs. continuous predictors is different.
For categorical predictors (Cat5), you can use emmeans() to give you the effects on the response scale directly:
emmOut <-emmeans(object = bestMod4.p, # your modelspecs =~ Cat5 + Cont6 + Cat5:Cont6 +1, # your predictorstype ="response") # report coefficients on the response scaleemmOut
Cat5 Cont6 rate SE df asymp.LCL asymp.UCL
StA 73.5 107.1 2.32 Inf 102.7 112
StB 73.5 95.2 2.61 Inf 90.2 100
StC 73.5 109.7 2.70 Inf 104.5 115
Confidence level used: 0.95
Intervals are back-transformed from the log scale
When Cat5 is StA, Resp4.p is estimated to be 74 (95% confidence interval: to 103 units).
When Cat5 is StB, Resp4.p is estimated to be 74 (95% confidence interval: to 90 units). When Cat5 is StC, Resp4.p is estimated to be 74 (95% confidence interval: to 105 units).
For continuous predictors (Cont6), you need to use the emtrends() function and then convert the coefficients to the response scale:
trendsOut <-emtrends(bestMod4.p,specs =~ Cat5 + Cont6 + Cat5:Cont6 +1, # your predictorsvar ="Cont6") # your predictor of interesttrendsOut
Since you have a poisson error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with the exp() function:
trendsOut.df <-data.frame(trendsOut) # convert trendsOut to a dataframeRR <-exp(trendsOut.df$Cont6.trend) # get the coefficient on the response scaleRR
[1] 1.004084 1.001458 1.008005
You can also use the same conversion on the confidence limits of the modelled effect, for the upper:
RR.up <-exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scaleRR.up
[1] 1.005716 1.004309 1.009791
and lower 95% confidence level:
RR.low <-exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scaleRR.low
[1] 1.0024542 0.9986158 1.0062224
Let’s organize these numbers so we can read the effects easier:
RRAll <-data.frame(Cat5 = trendsOut.df$Cat5, # include info on the Cat5 levelRR =exp(trendsOut.df$Cont6.trend), # the modelled effect as a rate ratioRR.up =exp(trendsOut.df$asymp.UCL), # upper 95% confidence level of the modelled effectRR.down =exp(trendsOut.df$asymp.LCL)) # lower 95% confidence level of the modelled effectRRAll
This tells you that for a unit change in Cont6, Resp4.p
increases by 1.004 times (an increase of 0.41%) when Cat5 is StA (95% confidence interval: 1.002 to 1.006).
increases by 1.001 times (an increase of 0.15%) when Cat5 is StB (95% confidence interval: 0.999 to 1.004).
increases by 1.008 times (an increase of 0.8%) when Cat5 is StC (95% confidence interval: 1.006 to 1.01).
Note that the 95% confidence interval for the effect of Cont6 on Resp4.p (the rate ratio) includes 1. When the rate ratio is 1, there is no effect of the continuous predictor on the response. See Section 5.2.4 for more explanation on this.
You can test for evidence that the effects of Cont6 on Resp4.p are significant with:
test(trendsOut) # get test if coefficient is different than zero.
Note that these tests are done on the link scale but can be used for reporting effects on the response scale.
From the results, you can see that:
there is evidence that there is an effect of Cont6 on Resp4.p when Cat5 = StA. (i.e. the slope on the link scale is different than zero, and the rate ratio on the response scale is different than 1; P < 0.0001).
there is no evidence that there is an effect of Cont6 on Resp4.p when Cat5 = StB. (i.e. the slope on the link scale is not different than zero, and the rate ratio on the response scale is not different than 1; P = 0.32).
there is evidence that there is an effect of Cont6 on Resp4.p when Cat5 = StC. (i.e. the slope on the link scale is different than zero, and the rate ratio on the response scale is different than 1; P < 0.0001).
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
You can also find out which effects of Cont6 on Resp4.p are different from one another using pairs() on the emtrends() output:
Cont6 = 73.5:
contrast estimate SE df z.ratio p.value
StA - StB 0.00262 0.00167 Inf 1.568 0.2598
StA - StC -0.00390 0.00123 Inf -3.180 0.0042
StB - StC -0.00652 0.00171 Inf -3.814 0.0004
P value adjustment: tukey method for comparing a family of 3 estimates
Remember that you need to convert any predictor coefficients to the response scale when reporting.
The table above tells you that:
There is no evidence the effect of Cont6 on Resp4.p differs when Cat5 is StA vs. StB (P = 0.26).
There is evidence that the effect of Cont6 on Resp4.p is different when Cat5 is StA vs. StC (P = 0.0042). The effect of Cont6 on Resp4.p is exp(-0.0039) = 0.9961076 times as big when Cat5 is StA vs. StC.
There is evidence that the effect of Cont6 on Resp4.p is different when Cat5 is StB vs. StC (P = 0.0004). The effect of Cont6 on Resp4.p is exp(-0.00652) = 0.9935012 times as big when Cat5 is StB vs. StC.
Only if all the slopes were similar, you would want to test if the levels of your categorical predictor (Cat5) have significantly different coefficients (intercepts) from each other with pairs() on the output from emmeans() (the emmOut object).
5.2.5 Models with a Gamma error distribution assumption:
Background
For models with a Gamma error distribution assumption, you will typically start with a link function that is the inverse or link = "inverse". You need to take this link function into consideration when you report your modelled effects.
For a model where link = "inverse", the interpretation of the coefficients on the response scale is less clear. In this case, you would use your model to make predictions of your response and report the effects by describing these changes in predictions. We will cover making predictions in the next chapter.
However, you can also use link = "log" for your models with a Gamma error distribution assumption. Using “log” will let you use the reporting methods in the section on “Models with a poisson error distribution assumption”(Section 5.2.4). So, if you have a model with a Gamma error distribution assumption, use link = "log" for modelling, and follow the methods given above (Section 5.2.4).
5.2.6 Models with a binomial error distribution assumption:
Background
For models with a binomial error distribution assumption, you will typically start with a link function that is the logit function or link = "logit". You need to take this link function into consideration when you report your modelled effects.
The logit function is also called the “log-odds” function as it is equal to the logarithm of the odds.
Models with a binomial error distribution assumption model the probability of an event happening. This might be presence (vs. absence), mature (vs. immature), death (vs. alive), etc., but in all cases this is represented as the probability of your response being 1 (vs. 0).
If \(p\) is the probability of an event happening (e.g. response being 1), the probability an event doesn’t happen is \(1-p\) when an event can only have two outcomes. The odds are the probability of the event happening divided by the probability it doesn’t happen, or:
\(\frac{p}{1-p}\)
The logit function (or log-odds function) is then:
\(log_e(\frac{p}{1-p})\)
Jargon alert! Note that a GLM using the logit link function is also called logistic regression.
Alternate link functions for binomially distributed data are probit and cloglog, but the logit link function is by far the one most used as the coefficients that result from the model are the easiest to interpret. But, as always, you should choose the link function that leads to the best-behaved residuals for your model.
Similar to reporting a rate ratio when your link = “log” (Section 5.2.4), you report your coefficients of a model with a binomial error distribution assumption (and link = “logit”) on the response scale by raising e to the power of the coefficients estimated on the link scale.
When you have a continuous predictor, \(e\) is raised to the slope (estimated on the link scale) to give you an odds ratio (OR) on the response scale. 1 - OR gives you the % change in odds for a unit change in your continuous predictor.
When you have a categorical predictor, \(e\) is raised to the intercept of each level of your predictor (estimated on the link scale) to give you the odds associated with that predictor level on the response scale.
Here is the math linking the coefficient to the odds ratio:
the math
Given a binomial model, \[
\begin{align}
log_e(\frac{p_i}{1-p_i}) &= \beta_1 \cdot Cov_i + \beta_0 \\
Resp_i &\sim binomial(p_i)
\end{align}
\]\(p_i\) is the probability of success and \(1-p_i\) is the probability of failure, and the odds are \(\frac{p_i}{1-p_i}\).
Then the odds ratio (OR), or % change in odds due to a change in predictor is:
The emmeans() function will convert the coefficients (intercepts) from the link to the response scale. You can ask for this with the type = "response" argument.
In contrast, the emtrends() function does not convert the coefficients (slopes) to represent effects on the response scale. This is because emtrends() is reporting the slope of a straight line - the trend line on the link scale. But the line isn’t straight on the response scale.
You need to convert from the link to the response by hand when we use emtrends()23. How you make the conversion, and interpret the result, will depend on your error distribution assumption:
Here we will go through examples of reporting for models with a binomial error distribution assumption. If you want more details on why you are using the methods below, check the section on “Models with a normal error distribution assumption” (Section 5.2.3) for context.
Examples
Here I have made new examples where the error distribution assumption is binomial This changes the response variables, indicated by Resp1.b, Resp2.b, etc. to be values of either 0 or 1.
Example 1: Resp1.p ~ Cat1 + 1
Example 1: Resp1.b ~ Cat1 + 1
Recall that Example 1 contains only one predictor and it is categorical:
Resp1.b ~ Cat1 + 1
The best-specified model (now with a binomial error distribution assumption) is:
bestMod1.b
Call: glm(formula = Resp1.b ~ Cat1 + 1, family = binomial(link = "logit"),
data = myDat1)
Coefficients:
(Intercept) Cat1StB Cat1StC
-0.6931 -1.8718 1.6094
Degrees of Freedom: 49 Total (i.e. Null); 47 Residual
Null Deviance: 68.03
Residual Deviance: 51.43 AIC: 57.43
and the dredge() table you used to pick your bestMod1.b is in dredgeOut1.b
dredgeOut1.b
Global model call: glm(formula = Resp1.b ~ Cat1 + 1, family = binomial(link = "logit"),
data = myDat1)
---
Model selection table
(Intrc) Cat1 R^2 df logLik AICc delta weight
2 -0.6931 + 0.2825 3 -25.714 57.9 0.00 0.998
1 -0.3228 0.0000 1 -34.015 70.1 12.16 0.002
Models ranked by AICc(x)
And here’s a visualization of the model effects:
# Set up your predictors for the visualized fitforCat1<-unique(myDat1.b$Cat1) # every value of your categorical predictor# create a data frame with your predictorsforVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of predictors# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod1.b, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data framelibrary(ggplot2)ggplot() +# start ggplotgeom_point(data = myDat1.b, # add observations to your plotmapping =aes(x = Cat1, y = Resp1.b)) +geom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),linewidth=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat1, y = Fit), shape =22, # shape of pointsize =3, # size of pointfill ="white", # fill colour for plotcol ='black') +# outline colour for plotylab("probability Resp1.b = 1")+# change y-axis labelxlab("Cat1, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
For categorical predictors, you can get the coefficients on the response scale directly with the emmeans() function24 when you specify type = “response”. These are the probabilities (not odds or odds ratio) that Resp1.b = 1 at each level of your categorical predictor:
library(emmeans) # load the emmeans packageemmOut <-emmeans(object = bestMod1.b, # your modelspecs =~ Cat1, # your predictorstype ="response") # report coefficients on the response scaleemmOut
Cat1 prob SE df asymp.LCL asymp.UCL
StA 0.3333 0.1220 Inf 0.14596 0.594
StB 0.0714 0.0688 Inf 0.00996 0.370
StC 0.7143 0.0986 Inf 0.49239 0.866
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
Notice that the column reported is prob: these coefficients are already converted back onto the response scale (as probabilities).
Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if you observed a new Resp1.b at a particular Cat1, there is a 95% probability that the probability Resp1.b is 1 would lie within the bounds of the confidence limits.
From this output you see:
When Cat1 is SpA, the probability that Resp1.b is 1 is estimated to be 0.33 (95% confidence interval: 0.15 to 0.59 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 0.5.
When Cat1 is SpB, the probability that Resp1.b is 1 is estimated to be 0.07 (95% confidence interval: 0.01 to 0.37 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 0.08.
When Cat1 is SpC, the probability that Resp1.b is 1 is estimated to be 0.71 (95% confidence interval: 0.49 to 0.87 units). This is equivalent to odds (\(\frac{p}{(1-p)}\)) of 2.5.
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
where do these probabilities and odds come from?
The emmeans package provides an easy way of seeing the effects of Cat1 on Resp1.b the response scale, but you can also get this information from the summary() output:
You can convert the coefficients to odds with the exp() function, for example:
When Cat1 is SpA, the odds that Resp1.b is 1 is estimated to be exp(coef(summary(bestMod1.b))[1]) = 0.5. And you can get the probabilities directly with the plogis() function (plogis(coef(summary(bestMod1.b))[1]) = 0.333).
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
Here you can compare the effects of the levels of Cat1 on Resp1.b to see if the effects differ between the levels.
forComp <-pairs(emmOut)forComp
contrast odds.ratio SE df null z.ratio p.value
StA / StB 6.5000 7.6300 Inf 1 1.595 0.2477
StA / StC 0.2000 0.1460 Inf 1 -2.204 0.0705
StB / StC 0.0308 0.0352 Inf 1 -3.041 0.0067
P value adjustment: tukey method for comparing a family of 3 estimates
Tests are performed on the log odds ratio scale
The output shows the results of the multiple comparison (pairwise) testing.
The values in the p.value column tell you the results of the hypothesis test comparing the coefficients between the two levels of your categorical predictor.
Note the odds.ratio column here. Comparing effects of Cat1 on the probability of Resp1.b is done by comparing odds of Resp1.b being 1 under each level of Cat1.
For example:
there is no evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpA vs. when Cat1 = SpB (P = 0.2477).
there is little evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpA than when Cat1 = SpC (P = 0.0705). The odds ratio is 0.2.
there is strong evidence that the probability of Resp1.b being 1 is different when Cat 1 = SpB than when Cat1 = SpC (P = 0.0067). The odds ratio is 0.0308.
where did this odds ratio come from?
This is a ratio of the odds of that Resp1.b is 1 when Cat 1 = SpA vs. the odds of that Resp1.b is 1 when Cat 1 = SpC:
For example, the odds that Resp1.b is 1 when Cat 1 = SpA are:
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCat2<-unique(myDat2.b$Cat2) # every level of your categorical predictorforCat3<-unique(myDat2.b$Cat3) # every level of your categorical predictor# create a data frame with your predictorsforVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of predictors#### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod2.b, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat2.b, # add observations to your plotmapping =aes(x = Cat2, y = Resp2.b, col = Cat3), position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) +# control position of data points so they are easier to see on the plotgeom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plotsize=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat2, y = Fit, fill = Cat3), shape =22, # shape of pointsize=3, # size of pointcol ='black', # outline colour for pointposition=position_dodge(width=0.75)) +# control position of data points so they are easier to see on the plotylab("probability Resp2.b = 1")+# change y-axis labelxlab("Cat2, (units)")+# change x-axis labellabs(fill="Cat3, (units)", col="Cat3, (units)") +# change legend titletheme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
As above, we can use the emmeans package to see the effects of the predictors on the scale of the response:
emmOut <-emmeans(object = bestMod2.b, # your modelspecs =~ Cat2 + Cat3 + Cat2:Cat3, # your predictorstype ="response") # report coefficients on the response scaleemmOut
Cat2 Cat3 prob SE df asymp.LCL asymp.UCL
StA Control 0.462 0.0523 Inf 0.362 0.564
StB Control 0.521 0.0585 Inf 0.407 0.632
StC Control 0.441 0.0515 Inf 0.344 0.543
StA Treat 0.436 0.0561 Inf 0.331 0.547
StB Treat 0.422 0.0521 Inf 0.325 0.526
StC Treat 0.840 0.0423 Inf 0.739 0.907
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
So now you have a modelled value of your response for each level of your categorical predictors. For example:
When Cat2 is StA and Cat3 is Treat, the probability that Resp1.b is 1 is estimated to be 0.44 (95% confidence interval: 0.33 to 0.55 units).
When Cat2 is StB and Cat3 is Control, the probability that Resp1.b is 1 is estimated to be 0.52 (95% confidence interval: 0.41 to 0.63 units).
Finally, note that you can also get a quick plot of the effects by handing the emmeans() output to plot().
Are the modelled effects different than one another across categorical levels? (only needed for categorical predictors with more than two levels)
As with Example 1, you can also find out which combinations of predictor levels are leading to statistically different model predictions in Resp2.b:
forComp <-pairs(emmOut, # emmeans outputsimple ="Cat2") # contrasts within this categorical predictorforComp
Cat3 = Control:
contrast odds.ratio SE df null z.ratio p.value
StA / StB 0.789 0.2490 Inf 1 -0.751 0.7331
StA / StC 1.087 0.3220 Inf 1 0.282 0.9572
StB / StC 1.377 0.4320 Inf 1 1.019 0.5646
Cat3 = Treat:
contrast odds.ratio SE df null z.ratio p.value
StA / StB 1.057 0.3300 Inf 1 0.179 0.9826
StA / StC 0.147 0.0573 Inf 1 -4.925 <.0001
StB / StC 0.139 0.0530 Inf 1 -5.183 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates
Tests are performed on the log odds ratio scale
Based on a threshold P-value of 0.05, you can see that the effects at some combinations of your predictors are not statistically different from each other. For example, there is no evidence that the probability of Resp2.b being 1 when Cat3 = Control and Cat2 = StA (0.46)25 than when Cat3 = Control and Cat2 = StB (0.52) are different from each other (P = 0.73).
On the other hand, some combinations of your predictors are statistically different from each other. For example, there is strong evidence that the probability of Resp2.b being 1 when Cat3 = Treat and Cat2 = StA (0.46) is less than that when Cat3 = Treat and Cat2 = StC (0.84) are different from each other (P < 0.0001).
Note that you can also get the results from the pairwise testing visually by handing the output of pairs() to plot().
Example 3: Resp3.p ~ Cont4 + 1
Example 3: Resp3.b ~ Cont4 + 1
Recall that Example 3 contains one predictor and the predictor is continuous:
Resp3 ~ Cont4 + 1
The best-specified model (now with a binomial error distribution assumption) is:
bestMod3.b
Call: glm(formula = Resp3.b ~ Cont4 + 1, family = binomial(link = "logit"),
data = myDat3.b)
Coefficients:
(Intercept) Cont4
-8.1519 0.1086
Degrees of Freedom: 299 Total (i.e. Null); 298 Residual
Null Deviance: 415.8
Residual Deviance: 199.1 AIC: 203.1
that was fit to data in myDat3.b:
str(myDat3.b)
'data.frame': 300 obs. of 2 variables:
$ Cont4 : num 57.7 53.2 79.7 35.1 72.2 ...
$ Resp3.b: int 1 0 1 0 0 1 1 1 1 0 ...
and the dredge() table you used to pick your bestMod3.b is
dredgeOut3.b
Global model call: glm(formula = Resp3.b ~ Cont4 + 1, family = binomial(link = "logit"),
data = myDat3.b)
---
Model selection table
(Intrc) Cont4 R^2 df logLik AICc delta weight
2 -8.15200 0.1086 0.5144 2 -99.533 203.1 0.00 1
1 0.04001 0.0000 1 -207.884 417.8 214.68 0
Models ranked by AICc(x)
And here’s a visualization of the model effects:
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCont4<-seq(from =min(myDat3.b$Cont4), to =max(myDat3.b$Cont4), by =1)# a sequence of numbers in your continuous predictor range# create a data frame with your predictorsforVis<-expand.grid(Cont4 = forCont4) # expand.grid() function makes sure you have all combinations of predictors. #### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod3.b, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat3.b, # add observations to your plotmapping =aes(x = Cont4, y = Resp3.b)) +# control position of data points so they are easier to see on the plotgeom_line(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cont4, y = Fit),size =1.2) +# control thickness of linegeom_line(data = forVis, # add uncertainty to your plot (upper line)mapping =aes(x = Cont4, y = Upper),size =0.8, # control thickness of linelinetype =2) +# control style of linegeom_line(data = forVis, # add uncertainty to your plot (lower line)mapping =aes(x = Cont4, y = Lower),size =0.8, # control thickness of linelinetype =2) +# control style of lineylab("probability Resp2.b = 1")+# change y-axis labelxlab("Cont4, (units)") +# change x-axis labeltheme_bw() # change ggplot theme
The model fit line may seem a little strange as the data are binary (can only be 0 or 1) but the model fit line goes through the space in between 0 and 1. This model fit line tells you how the probability of Resp3.b being 1 changes as Cont4 changes.
What are your modelled effects (with uncertainty)?
trendsOut <-emtrends(bestMod3.b,specs =~ Cont4 +1, # your predictorsvar ="Cont4") # your predictor of interesttrendsOut
Note that this is the same estimate as you find in the summary() output of your model object:
coef(summary(bestMod3.b))
Estimate Std. Error z value Pr(>|z|)
(Intercept) -8.1518996 0.89160247 -9.142976 6.075711e-20
Cont4 0.1085946 0.01160302 9.359159 8.037377e-21
This is because emtrends() returns coefficients on the link scale.
To report the coefficients on a response scale, you need to consider your link function. Since you have a binomial error distribution assumption, you can convert the estimate made by emtrends() to the response scale with:
trendsOut.df <-data.frame(trendsOut) # convert trendsOut to a dataframeOR <-exp(trendsOut.df$Cont4.trend) # get the coefficient on the response scaleOR
[1] 1.11471
You can also use the same conversion on the confidence limits of the modelled effect, for the upper:
OR.up <-exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scaleOR.up
[1] 1.140351
and lower 95% confidence level:
OR.low <-exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scaleOR.low
[1] 1.089646
This coefficient - called the odds ratio (Section 5.2.6) - tells you that for a unit change in Cont4, the odds that Resp3.b is 1 increases by 1.115 times (95% confidence interval is 1.09 to 1.14) which is an increase in the odds of Resp3.b being 1 of 11.5%.
where did this odds ratio come from?
You can get a better sense of the odds ratio (vs. odds vs. probabilities) by estimating it directly from the modelled probability that Resp3.b is 1:
Let’s find the probability of Resp3.b = 1 when Cont4 is 60:
pCont4.60<-predict(bestMod3.b, # modelnewdata =data.frame(Cont4 =60), # predictor values for the predictiontype ="response", # give prediction on the response scalese.fit =TRUE) # include uncertaintypCont4.60<- pCont4.60$fitpCont4.60
1
0.1629792
So there is a 16.3% probability of Resp3.b being 1 when Cont4 is 60.
The odds associated with Resp3.b being 1 when Cont4 is 60 is given by the probability of Cont4 is divided by the probability of absence:
Now let’s find the probability of Resp3.b being 1 when Cont4 is 61:
pCont4.61<-predict(bestMod3.b, # modelnewdata =data.frame(Cont4 =61), # predictor values for the predictiontype ="response", # give prediction on the response scalese.fit =TRUE) # include uncertaintypCont4.61<- pCont4.61$fitpCont4.61
1
0.1783404
So there is a 17.8% probability of Resp3.b being 1 when Cont4 is 61.
The odds associated with Resp3.b being 1 when Cont4 is 61 is given by the probability of presence divided by the probability of absence:
#### i) choosing the values of your predictors at which to make a prediction# Set up your predictors for the visualized fitforCat5<-unique(myDat4.b$Cat5) # all levels of your categorical predictorforCont6<-seq(from =min(myDat4.b$Cont6), to =max(myDat4.b$Cont6), by =1)# a sequence of numbers in your continuous predictor range# create a data frame with your predictorsforVis<-expand.grid(Cat5 = forCat5, Cont6 = forCont6) # expand.grid() function makes sure you have all combinations of predictors. #### ii) using `predict()` to use your model to estimate your response variable at those values of your predictor# Get your model fit estimates at each value of your predictorsmodFit<-predict(bestMod4.b, # the modelnewdata = forVis, # the predictor valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame#### iii) use the model estimates to plot your model fitlibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat4.b, # add observations to your plotmapping =aes(x = Cont6, y = Resp4.b, col = Cat5)) +# control position of data points so they are easier to see on the plotgeom_line(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cont6, y = Fit, col = Cat5),size =1.2) +# control thickness of linegeom_line(data = forVis, # add uncertainty to your plot (upper line)mapping =aes(x = Cont6, y = Upper, col = Cat5),size =0.4, # control thickness of linelinetype =2) +# control style of linegeom_line(data = forVis, # add uncertainty to your plot (lower line)mapping =aes(x = Cont6, y = Lower, col = Cat5),size =0.4, # control thickness of linelinetype =2) +# control style of lineylab("probability Resp2.b = 1")+# change y-axis labelxlab("Cont6, (units)") +# change x-axis labellabs(fill="Cat5, (units)", col="Cat5, (units)") +# change legend titletheme_bw() # change ggplot theme
What are your modelled effects (with uncertainty)?
The process for estimating effects of categorical vs. continuous predictors is different.
For categorical predictors (Cat5), you can use emmeans() to give you the effects on the response scale directly:
emmOut <-emmeans(object = bestMod4.b, # your modelspecs =~ Cat5 + Cont6 + Cat5:Cont6 +1, # your predictorstype ="response") # report coefficients on the response scaleemmOut
Cat5 Cont6 prob SE df asymp.LCL asymp.UCL
StA 198 0.951 0.0346 Inf 0.819 0.988
StB 198 0.964 0.0178 Inf 0.907 0.986
StC 198 0.438 0.0941 Inf 0.269 0.622
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
This output tells you that:
When Cat5 is StA and Cont6 is 227 units26, there is a 95% probability that Resp4.b will be 1 (95% confidence interval: 0.82 to 0.99%)
When Cat5 is StB and Cont6 is 227 units27, there is a 96% probability that Resp4.b will be 1 (95% confidence interval: 0.91 to 0.99 units).
When Cat5 is StB and Cont6 is 227 units28, there is a 44% probability that Resp4.b will be 1 (95% confidence interval: 0.27 to 0.62 units).
Note that emmeans() sets your continuous predictor (Cont6) to the mean value of Cont6 (227 units). You can also control this with the at = argument in the emmeans() function:
emmOut <-emmeans(object = bestMod4.b, # your modelspecs =~ Cat5 + Cont6 + Cat5:Cont6 +1, # your predictorstype ="response", # report coefficients on the response scaleat =list(Cont6 =150)) # control the value of your continuous predictor at which to make the coefficient estimatesemmOut
Cat5 Cont6 prob SE df asymp.LCL asymp.UCL
StA 150 8.93e-03 0.009230 Inf 1.17e-03 0.06500
StB 150 5.53e-01 0.078500 Inf 3.99e-01 0.69695
StC 150 6.27e-05 0.000142 Inf 7.00e-07 0.00524
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
For continuous predictors (Cont6), you need to use the emtrends() function and then convert the coefficients to the response scale:
trendsOut <-emtrends(bestMod4.b,specs =~ Cat5 + Cont6 + Cat5:Cont6 +1 , # your predictorsvar ="Cont6") # your predictor of interesttrendsOut
Since you have a binomial error distribution assumption, you can convert the estimate made by emtrends() on to the response scale with the exp() function:
trendsOut.df <-data.frame(trendsOut) # convert trendsOut to a dataframeOR <-exp(trendsOut.df$Cont6.trend) # get the coefficient on the response scaleOR
[1] 1.172328 1.065611 1.215576
You can also use the same conversion on the confidence limits of the modelled effect, for the upper:
OR.up <-exp(trendsOut.df$asymp.UCL) # get the upper confidence interval on the response scaleOR.up
[1] 1.250737 1.089631 1.325976
and lower 95% confidence level:
OR.low <-exp(trendsOut.df$asymp.LCL) # get the lower confidence interval on the response scaleOR.low
[1] 1.098834 1.042120 1.114368
Let’s organize these numbers so we can read the effects easier:
ORAll <-data.frame(Cat5 = trendsOut.df$Cat5, # include info on the Cat5 levelOR =exp(trendsOut.df$Cont6.trend), # the modelled effect as a rate ratioOR.up =exp(trendsOut.df$asymp.UCL), # upper 95% confidence level of the modelled effectOR.down =exp(trendsOut.df$asymp.LCL)) # lower 95% confidence level of the modelled effectORAll
Cat5 OR OR.up OR.down
1 StA 1.172328 1.250737 1.098834
2 StB 1.065611 1.089631 1.042120
3 StC 1.215576 1.325976 1.114368
This tells you that for a unit change in Cont6, the odds that Resp4.b is 1:
increases by 17.2% with a unit change in Cont6 when Cat5 is StA (95% confidence interval: 9.9 to 25.1).
increases by 6.6% with a unit change in Cont6 when Cat5 is StB (95% confidence interval: 4.2 to 9).
increases by 21.6% with a unit change in Cont6 when Cat5 is StC (95% confidence interval: 11.4 to 32.6).
You can test for evidence that the effects of Cont6 on Resp4.b are significant (different than zero) with:
test(trendsOut) # get test if coefficient is different than zero.
Cont6 = 198:
contrast estimate SE df z.ratio p.value
StA - StB 0.0954 0.0349 Inf 2.732 0.0173
StA - StC -0.0362 0.0553 Inf -0.655 0.7895
StB - StC -0.1317 0.0458 Inf -2.876 0.0112
P value adjustment: tukey method for comparing a family of 3 estimates
The table above tells you that:
There is evidence that the effect of Cont6 on Resp4.b is different when Cat5 is StA vs. StB (P = 0.017). The effect of Cont6 on Resp4.b is exp(0.0954) = 1.1 times as big when Cat5 is StA vs. StB (i.e. the effect of Cont6 on Resp4.b is 10% more when Cat5 is StA vs. StB)
There is no evidence the effect of Cont6 on Resp4.b differs when Cat5 is StA vs. StC (P = 0.789).
There is evidence that the effect of Cont6 on Resp4.b is different when Cat5 is StB vs. StC (P = 0.011). The effect of Cont6 on Resp4.b is exp(-0.1317) = 0.88 times as big when Cat5 is StB vs. StC (i.e. the effect of Cont6 on Resp4.b is 12% less when Cat5 is StB vs. StC)
Only if all the slopes29 were similar, you would want to test if the levels of your categorical predictor (Cat5) have significantly different coefficients (intercepts) from each other with pairs() on the output from emmeans() (the emmOut object).
Up next
In the next chapter, we will discuss how you can use your model to make predictions.
6 Communicating your hypothesis testing results
Here we’ll walk through your statistical modelling framework and discuss how to communicate motivations, methods and results of your data analysis and hypothesis testing for a biological report/article (& the exam!).
Here again is our statistical modelling framework:
Note that I’ve merged the Visualizing section with the Reporting section to correspond with how we will write-up our results.
6.1 A familiar example
I’ll use the following example case as I discuss how to communicate your results:
You are interested in moose (or elk, Alces alces) and how their growth rate changes from population to population. You find information on mean annual weight change by moose in 12 locations in Sweden. You wonder if latitude can explain some of the variability in growth rates you are observing. You also wonder if minimum winter temperature might explain some variability in growth rate. Finally, you wonder if changes in growth rate with latitude and/or winter temperature depend on each other (i.e, latitudinal dependent effects on growth rate depend on what the minimum winter temperature is) and/or sex (male and female) - i.e. the environmentally dependent change in growth rate depends on what sex the animal is.
And here’s a screenshot of the readme file for the moose data:
6.2 Response(s)
The task:
Your response variable30 is the observed variability you are trying to explain. Your response forms your research question - i.e. “why is my response varying?”. In this section, you need to identify the variability you’re trying to explain and why (your motivation).
The code:
Code supporting your tasks in this section allow you to import your data and help you describe your response variable, e.g.:
rm(list=ls()) # make sure your workspace is clear#Import the datamyDat <-read.table("ExDat.csv", # file namesep =',', # lines in the file are separated by a comma - this information is in the README filedec ='.', # decimals are denoted by '.' - this information is in the README filestringsAsFactors =TRUE, # make sure strings (characters) are imported as factorsheader =TRUE) # there is a column header in the filestr(myDat) # check the structure of the data file
summary(myDat) # get summary information of each column
Study Sex delMass Lat WTemp
Min. : 1.00 Females:12 Min. :36.20 Min. :57.70 Min. :-11.900
1st Qu.: 3.75 Males :12 1st Qu.:41.95 1st Qu.:58.00 1st Qu.: -8.400
Median : 6.50 Median :47.35 Median :61.80 Median : -6.600
Mean : 7.00 Mean :47.77 Mean :61.60 Mean : -6.100
3rd Qu.: 9.75 3rd Qu.:54.23 3rd Qu.:64.38 3rd Qu.: -2.225
Max. :14.00 Max. :60.60 Max. :66.00 Max. : -1.700
Source
Gothenburg: 6
Lund : 8
Stockholm :10
The write-up:
Use the code above, the data readme file, as well as supporting literature to communicate about your response variable and motivation. Include:
what your response variable is
why you want to explain variability in your response
your research question
how your response variable was measured and a description of how your response varies. For example, give the range in your response if possible. If your response is a category, give the levels of the category (e.g. alive vs. dead). Always give units, if applicable.
For example:
Here I want to explain variability in growth rate changes in moose (Alces alces). Moose size will be related to survival rates, and likely also their ability to reproduce (e.g. Sand et al. 1995). Being able to explain what controls growth rate in moose may help us explain why some populations are more productive than other populations, and it will help us predict productivity for populations at other places or times.
My research question is “why does growth rate vary?”. My observations of growth rate were measured from adult moose across 12 locations in Sweden. The growth rate varied from 36.2 to 60.6 kg per year.
6.3 Covariate(s)
The task:
Here you will identify the covariates in your model - i.e. what variables might explain variability in your response. To write-up this section, start with the biological mechanism (or processes) that make you think the covariate could be influencing your response. Then describe how that covariate was measured. At this point you might identify mechanisms that you are not able to measure, e.g. because they were beyond the scope of your study. Make note of these. You will have a chance to talk about them at the end of this exercise (in Reporting).
The code:
Code supporting your tasks in this section will help you describe your covariates, e.g.:
# Getting information on my covariatesrange(myDat$Lat) # ranges for continuous covariates
[1] 57.7 66.0
range(myDat$WTemp) # ranges for continuous covariates
[1] -11.9 -1.7
unique(myDat$Sex) # factor levels for categorical covariates
[1] Females Males
Levels: Females Males
The write-up:
Use the code above, the data readme file, as well as supporting literature to communicate about your covariate variables. Include:
what your covariate variables are
why you think each covariate might explain variability in your response (mechanisms)
how each covariate variable was measured and a description of the covariate. This means at least giving the range of each covariate (for continuous covariates) or the levels of any categorical covariate. Remember to include units where applicable.
For example:
Here I consider the effects of latitude (˚N), minimum winter temperature (˚C) and sex (male vs. female) on annual growth in adult moose.
Latitude may impact the growth of the moose as the environment changes with latitude (e.g. temperature, food, light level, growing season). Differences in food availability may lead to growth differences [citation]31. In this study, latitude is measured in decimal degrees N and ranges from 57.7 to 66.0˚N.
Winter harshness may impact annual moose growth by making energetic losses larger in the winter due to increased costs of temperature regulation [citation]. Here, I measure winter harshness as minimum winter temperature (˚C) which ranged from -11.9 to -1.7˚C in my study.
Growth differences might also occur due to sex of the animal as male and female moose have different life history strategies (e.g. reproductive investment)[citation]. Here, I include effects of sex with observations of male vs. female for each growth rate measurement.
6.4 Hypothesis
The task:
In this section you will describe your research hypothesis including main effects and any interactions. Remember to include the formula notation we’ve been practicing in class (i.e. Response ~ Covariate…)
The code:
No code needed here!
The write-up:
Use your answers to the previous sections to state your research hypothesis. Remember to consider interactions as well as the main effects. Include:
a description in words
a description in the formula notation
a definition of each term
For example,
I will test the research hypothesis that varibility in adult moose growth rate (delMass, kg per year) is explained by latitude (Lat, ˚N), minimum winter temperature (WTemp, ˚C) and sex (Sex, male or female). I will test this by modelling:
delMass ~ Lat + WTemp + Sex + Lat:Sex + WTemp:Sex + Lat:WTemp + Lat:WTemp:Sex
which includes main effects and all possible interaction effects of my covariates (where “:” indicates an interaction between covariates).
6.5 Starting model
The task:
This step involves describing how you built the statistical model to test your hypothesis. Here, you will:
Choose your error distribution assumption
Choose your shape assumption
Choose your starting model
Fit your starting model
The code:
Code supporting your tasks in this section will help you support your choices of your model assumptions, e.g.:
# Choosing your error distribution assumption:## Understanding the nature of your response variablesummary(myDat$delMass) # a summary description of delMass
Min. 1st Qu. Median Mean 3rd Qu. Max.
36.20 41.95 47.35 47.77 54.23 60.60
# Choosing your shape assumption:## This can include plotting your response variable vs. each continuous covariate### For Lat:library(ggplot2) # load ggplot2 libraryggplot()+# start ggplotgeom_point(data = myDat,mapping =aes(x = Lat, y = delMass))+# add observationsxlab("Latitude, ˚N")+# change x-axis labelylab("Annual growth rate, kg per year")+# change y-axis labellabs(caption ="Figure 1: Observed annual growth rate (kg per year) vs. latitude (˚N)")+# add figure captiontheme_bw()+# change theme of plottheme(plot.caption =element_text(hjust=0)) # move figure legend (caption) to left alignment. Use hjust = 0.5 to align in the center.
#### Relationship between Lat and delMass looks linear### For WTemp:ggplot()+# start ggplotgeom_point(data = myDat,mapping =aes(x = WTemp, y = delMass))+# add observationsxlab("Minimum winter temperature, ˚C")+# change x-axis labelylab("Annual growth rate, kg per year")+# change y-axis labellabs(caption ="Figure 2: Observed annual growth rate (kg per year) vs. minimum winter temperature (˚C)")+# add figure captiontheme_bw()+# change theme of plottheme(plot.caption =element_text(hjust=0)) # move figure legend (caption) to left alignment. Use hjust = 0.5 to align in the center.
#### Relationship between WTemp and delMass looks linear# Fitting your starting modelstartMod<-glm(formula = delMass ~ Lat + WTemp + Sex + Lat:Sex + WTemp:Sex + Lat:WTemp + Lat:WTemp:Sex, # hypothesisdata = myDat, # datafamily =gaussian(link ="identity")) # error distribution assumptionstartMod
Call: glm(formula = delMass ~ Lat + WTemp + Sex + Lat:Sex + WTemp:Sex +
Lat:WTemp + Lat:WTemp:Sex, family = gaussian(link = "identity"),
data = myDat)
Coefficients:
(Intercept) Lat WTemp SexMales
13.86368 0.43887 1.37886 -54.84374
Lat:SexMales WTemp:SexMales Lat:WTemp Lat:WTemp:SexMales
1.09654 -5.46475 -0.02401 0.08799
Degrees of Freedom: 23 Total (i.e. Null); 16 Residual
Null Deviance: 1200
Residual Deviance: 75.09 AIC: 113.5
The write-up:
Use the code in the previous section as well as your theoretical understanding of the variables to communicate how you chose and fit your starting model. Include:
what your error distribution assumption is and how you chose it
what your shape assumption is and how you chose it
what type of model you chose for your starting model (e.g. GLM)
For example,
I built a model to test my hypothesis that varibility in adult moose growth rate (delMass, kg per year) is explained by latitude (Lat, ˚N), minimum winter temperature (WTemp, ˚C) and sex (Sex, male or female). My error distribution assumption was normal (gaussian) as delMass is continuous and could be positive or negative. My shape assumption regarding the relationship between delMass and each covariate was linear. This was chosen after inspecting plots of delMass and both WTemp and Lat (Figures 1 & 2). Plots were made using the ggplot2 package (Wickham 2016)32 In addition, a linear shape assumption was chosen to describe the relationship between Sex and delMass as Sex is a categorical variable.
Based on these assumption I tested my hypothesis by fitting a generalized linear model with normal error distribution assumption to my data. All model fitting and analysis was done in R (R Core Team 2022).33
6.6 Assessing fit
The task:
Here, you will describe how you assessed your model to ensure that it can be used to test your hypothesis. Remember that the assumptions you made above (error distribution & shape assumptions) were a “best guess” but it is only after the model is fit that we can confirm that it is useful. In this section you will:
Consider covariate collinearity
Consider observation independence
Consider your error distribution assumption
Consider your shape assumption
and report your results.
The code:
Code supporting your tasks in this section will help you assess your model for violations that may make it invalid for testing your hypothesis, e.g.:
# Consider covariate collinearity ## Fit a model without interactionsstartMod.noInt <-glm(formula = delMass ~ Lat + WTemp + Sex, # hypothesisdata = myDat, # datafamily =gaussian(link ="identity")) # error distribution assumptionstartMod.noInt
Call: glm(formula = delMass ~ Lat + WTemp + Sex, family = gaussian(link = "identity"),
data = myDat)
Coefficients:
(Intercept) Lat WTemp SexMales
-13.013 0.879 -0.105 12.000
Degrees of Freedom: 23 Total (i.e. Null); 20 Residual
Null Deviance: 1200
Residual Deviance: 99.43 AIC: 112.2
## Calculate Variance Inflation Factors (VIFs)#library(car) # load the car package#vif(startMod.noInt) # estimate VIFs## there are VIF values higher than the threshold value (in this case VIFs > 3). Let's remove WTemp and recalculate the VIFs:startMod.noInt.noWTemp <-glm(formula = delMass ~ Lat + Sex, # hypothesisdata = myDat, # datafamily =gaussian(link ="identity")) # error distribution assumptionstartMod.noInt.noWTemp
Call: glm(formula = delMass ~ Lat + Sex, family = gaussian(link = "identity"),
data = myDat)
Coefficients:
(Intercept) Lat SexMales
-19.5975 0.9963 12.0000
Degrees of Freedom: 23 Total (i.e. Null); 21 Residual
Null Deviance: 1200
Residual Deviance: 99.58 AIC: 110.3
## Recalculating the VIFs#vif(startMod.noInt.noWTemp) # estimate VIFs## there are no longer any problems with covariate collinearity, as all VIFs < 3## Refitting your starting model with the interactions back in:startMod <-glm(formula = delMass ~ Lat + Sex + Lat:Sex, # hypothesisdata = myDat, # datafamily =gaussian(link ="identity")) # error distribution assumption# Consider observation independence## You wonder if Source could be violating your observation independence assumption. Check to see with a plot of your residuals vs. Sourcelibrary(DHARMa) # load packagesimulationOutput <-simulateResiduals(fittedModel = startMod, n =250) # simulate data from our model n times and calculate residualsmyDat$Resid <- simulationOutput$scaledResiduals # add residuals to data frameggplot()+# start ggplotgeom_violin(data = myDat,mapping =aes(x = Source, y = Resid))+# add observations as a violingeom_point(data = myDat,mapping =aes(x = Source, y = Resid))+# add observations as pointsxlab("Source university")+# y-axis labelylab("Scaled residual")+# x-axis labellabs(caption ="Figure 3: A comparison of model residuals vs. university from which the data were sourced")+# figure captiontheme_bw()+# change theme of plottheme(plot.caption =element_text(hjust=0)) # move figure legend (caption) to left alignment. Use hjust = 0.5 to align in the center.
## The pattern is uniform. You assume no problem with observation dependence.# Consider your error distribution assumption by inspecting residuals plotted vs. the fitted valuesplot(simulationOutput, asFactor=FALSE) # compare simulated data to our observations
# Consider your shape assumption by inspecting residuals plotted vs. each covariateplot(simulationOutput, # compare simulated data to form=myDat$Sex, # our observationsasFactor=TRUE) # whether the variable plotted is a factor
plot(simulationOutput, # compare simulated data to form=myDat$Lat, # our observationsasFactor=FALSE) # whether the variable plotted is a factor
The write-up:
Use the code in the previous section to comment on your model’s validity in testing your research hypothesis. Include:
how you determined if there were problems with covariate collinearity and any actions you took if you detected a problem
how you determined if there were problems with observation dependence and any actions you took if you detected a problem
how you determined if your error distribution assumption was valid and any actions you took to address problems
how you determined if your shape assumption was valid and any actions you took to address problems
For example,
I tested if collinearity among my covariates was making my model coefficients uncertain by estimating variance inflation factors (VIFs) with the car package (Fox & Weisberg 2019)34. Initial VIFs for Lat and WTemp were > 23 indicating a high level of covariate collinearity. WTemp was removed from the model and the VIFs were re-estimated. VIFs in the new model were both 1, and it was concluded that there was no further issue with covariate collinearity. The new starting model fits the hypothesis:
delMass ~ Latitude + Sex + Latitude:Sex
I tested my assumption of observation independence by determining if the observations were grouped by Source (the source university for the data35). I estimated my model residuals using the DHARMA package (Hartig 2022)36 and plotted my residuals vs. Source. I concluded my data were not dependent on one another (grouped by) Source as the residuals were evenly distributed across the three source universities (see figure 3).
I assessed my error distribution assumption by inspected my residuals. Observed residuals were similar to those expected given my normal error distribution assumption. The Kolmogorov-Smirnov test comparing the observed to the expected distribution was not significant (P = 0.96). The dispersion and outlier tests were also not significant (P = 0.74 and P = 0.99 respectively). A plot of the residuals vs. fitted values showed a uniform cloud. From these results, I concluded that my error distribution assumption was appropriate.
I assessed my shape assumption by inspected my residuals vs. each covariate. My residuals were evenly distributed with Latitude, indicating a linear shape assumption for Latitude was appropriate. The linear shape assumption for Sex was necessary as Sex is a categorical variable. A plot of the residuals vs. Sex showed residuals were evenly distributed across the two sexes. From these results, I concluded that my linear shape assumptions were appropriate.
Given these results, I determined that the new starting model (delMass ~ Latitude + Sex + Latitude:Sex) can be used to test my hypothesis.
6.7 Hypothesis test
The task:
Here, you will use your model to test your hypothesis. For the purposes of this course, you will do this with the model selection method we practiced in class37. In this section you will:
fit and compare models representing all possible combinations of the covariates in your starting model
use the results to choose best-specified model(s)
report your results.
The code:
Code supporting your tasks in this section will help you test and rank your models, e.g.:
library(MuMIn) # load packageoptions(na.action ="na.fail") # needed for dredge() function to prevent illegal model comparisons(dredgeOut<-dredge(startMod, extra ="R^2")) # fit and compare a model set representing all possible predictor combinations
Global model call: glm(formula = delMass ~ Lat + Sex + Lat:Sex, family = gaussian(link = "identity"),
data = myDat)
---
Model selection table
(Int) Lat Sex Lat:Sex R^2 df logLik AICc delta weight
8 -1.596 0.7041 + + 0.9340 5 -48.387 110.1 0.00 0.756
4 -19.600 0.9963 + 0.9170 4 -51.130 112.4 2.26 0.244
3 41.780 + 0.7200 3 -65.726 138.7 28.54 0.000
2 -13.600 0.9963 0.1971 3 -78.366 163.9 53.82 0.000
1 47.780 0.0000 2 -80.999 166.6 56.46 0.000
Models ranked by AICc(x)
## Our best-specified model is delMass ~ Lat + Sex + Lat:Sex (model #8 in the dredge() output).bestMod<-(eval(attr(dredgeOut, "model.calls")$`8`)) # extract model #8 from dredge table# You can make a pretty table to use in your write-uplibrary(gt) # load gt packagemyTable <-gt(dredgeOut) # make a pretty tablemyTable <-fmt_number(myTable, # to format the numbers in my tablecolumns =everything(), # which columns to formatdecimals =2) # round to 2 decimal placesmyTable <-tab_caption(myTable, caption ="Table 1: Model selection table for hypothesis testing. Each row is a model fit to my data. Covariates included in each model are indicated by a number (for continuous covariates and intercept) or '+' (for categorical covariates). R^2 is the likelihood ratio R-squared, df indicates number of model parameters, logLik is the model Log-likelihood, AICc is the corrected Akaike Information Criteria, delta is the change in AICc between the model and that of the lowest AICc, and weight is the Akaike weights.")myTable
Table 1: Model selection table for hypothesis testing. Each row is a model fit to my data. Covariates included in each model are indicated by a number (for continuous covariates and intercept) or '+' (for categorical covariates). R^2 is the likelihood ratio R-squared, df indicates number of model parameters, logLik is the model Log-likelihood, AICc is the corrected Akaike Information Criteria, delta is the change in AICc between the model and that of the lowest AICc, and weight is the Akaike weights.
(Intercept)
Lat
Sex
Lat:Sex
R^2
df
logLik
AICc
delta
weight
−1.60
0.70
+
+
0.93
5.00
−48.39
110.11
0.00
0.76
−19.60
1.00
+
NA
0.92
4.00
−51.13
112.36
2.26
0.24
41.78
NA
+
NA
0.72
3.00
−65.73
138.65
28.54
0.00
−13.60
1.00
NA
NA
0.20
3.00
−78.37
163.93
53.82
0.00
47.78
NA
NA
NA
0.00
2.00
−81.00
166.57
56.46
0.00
## Could also be written as:# dredgeOut %>%# gt() %>% # make a pretty table# fmt_number(columns = everything(), # which columns to format# decimals = 2) %>% # round to 2 decimal places# tab_caption(caption = "Table 1: Model selection table for hypothesis testing. Each row is a model fit to my data. Covariates included in each model are indicated by a number (for continuous covariates and intercept) or '+' (for categorical covariates). R^2 is the likelihood ratio R-squared, df indicates number of model parameters, logLik is the model Log-likelihood, AICc is the corrected Akaike Information Criteria, delta is the change in AICc between the model and that of the lowest AICc, and weight is the Akaike weights.")
The write-up:
Use the code in the previous section to explain how you chose your best-specified model. Include:
the method you are using to test your hypothesis
how you fit and rank your candidate model set
how you made your decision regarding your best-specified model
For example,
I used model selection to test my hypothesis that varibility in adult moose growth rate is explained by latitude, sex and the interaction between the two. I used the dredge() function from the MuMIn package (Bartón 2022) to fit and rank models representing all possible covariate combinations. Models were ranked by corrected Akaike Information Criteria38 which [add a brief description of what AICc is].39
The model with the lowest AICc was chosen as the best-specified model, with models within ∆AICc = 2 of the lowest AICc model being considered equally best-specified.
The best-specified model was the full model:
delMass ~ Lat + Sex + Lat:Sex
with an Akaike weight (normalized relative likelihood) of 0.756. The next best model had an AIC of 2.26 more than the top model and an Akaike weight of 0.244.
6.8 Reporting
The task:
In this section, you will report your hypothesis testing results. Specifically,
Report your best-specified model(s)
Report your modelled effects: This includes both reporting your coefficients and plotting the effects. I’ve moved visualization here to be more consistent with how we report biological results. A few things to note here:
reporting the coefficients includes:
reporting your coefficient estimates with uncertainty: make sure you report appropriate significant digits and units.
reporting evidence that these coefficient estimates are different than zero: the fact that your covariate is in your best-specified model means that there is evidence of an effect of the covariate on the response (i.e. that the coefficient associated with the covariate is different than zero), BUT if you have a categorical variable with multiple levels, you need to test to see which coefficients for each level are different than zero.
reporting evidence that these coefficient estimates are different than each other across the levels of your categorical covariate (if you have one): this is the multiple comparison testing we practiced in class. Remember that if you have both a continuous covariate and a categorical covariate, you need to test the coefficients (slopes) for the continuous covariate first and only if the slopes are similar across levels of the categorical covariate can you test if the intercepts are significantly different. This is due to the fact that we can get different intercept estimates just due to slopes not being equal.
In your visualization, remember to include your modelled fit, uncertainty and your observations in your plot)
Report how well your model explains your response (explained deviance)
Report the relative importance of each term in your model
The code:
Code supporting your tasks in this section will help you communicate what your results tell you about your hypothesis, e.g.:
# Reporting the modelled effects of your covariates on your response: ## ** your coefficient estimates with uncertainty### Find the coefficients for Latitude. As Sex is categorical, we do this for emmeans() and the reported coefficients are predicted mean response at each level of Sex. Note that emmeans() can report coefficients on the RESPONSE scale.library(emmeans) # load emmeans packageemmeans(object = bestMod, # your modelspecs =~ Lat + Sex + Lat:Sex, # your covariates# at = list(Lat = 0), # to get coefficient predictions for Sex when Lat = 0type ="response") # report coefficients on the response scale
Lat Sex emmean SE df lower.CL upper.CL
61.6 Females 41.8 0.575 20 40.6 43
61.6 Males 53.8 0.575 20 52.6 55
Confidence level used: 0.95
### Find the coefficients for Latitude. As Latitude is continuous, we do this with emtrends() and the reported coefficients are the slope describing the change in predicted mean response with a unit change in Latitude. Note that emtrends() reports coefficients on the LINK scale. You need to convert this to the response scale. In our case, the error distribution assumption is normal so nothing needs to be done to convert the coefficients.emtrends(object = bestMod, # your modelspecs =~ Lat + Sex + Lat:Sex, # your covariatesvar ="Lat") # name of your continuous cvoariate
Lat Sex Lat.trend SE df lower.CL upper.CL
61.6 Females 0.704 0.182 20 0.324 1.08
61.6 Males 1.289 0.182 20 0.908 1.67
Confidence level used: 0.95
## ** evidence that these coefficient estimates are different than zero & ** evidence that these coefficient estimates are different than each other across the levels of your categorical covariate (if you have one)### Are the coefficients (slopes) of the latitude effect on growth different?forLatCompare <-emtrends(object = bestMod, # your modelspecs = pairwise ~ Lat + Sex + Lat:Sex, # your covariatesvar ="Lat") # name of your continuous cvoariatetest(forLatCompare) # to get results of slope tests and comparisons
$emtrends
Lat Sex Lat.trend SE df t.ratio p.value
61.6 Females 0.704 0.182 20 3.861 0.0010
61.6 Males 1.289 0.182 20 7.065 <.0001
$contrasts
contrast estimate SE df t.ratio p.value
Lat61.6 Females - Lat61.6 Males -0.584 0.258 20 -2.266 0.0347
plot(forLatCompare, # plot coefficients with errorcomparison =TRUE, # include arrows indicating thresholds for differencestype ="response") # plot on the response scale
### Are the coefficients (intercepts) of Sex factor levels different than zero?forSexCompare <-emmeans(object = bestMod, # your modelspecs = pairwise ~ Lat + Sex + Lat:Sex, # your covariatestype ="response") # name of your continuous cvoariatetest(forSexCompare)
$emmeans
Lat Sex emmean SE df t.ratio p.value
61.6 Females 41.8 0.575 20 72.704 <.0001
61.6 Males 53.8 0.575 20 93.588 <.0001
$contrasts
contrast estimate SE df t.ratio p.value
Lat61.6 Females - Lat61.6 Males -12 0.813 20 -14.768 <.0001
## Plot your modelled effects:# Set up your covariates for the visualized fitforLat<-seq(from =min(myDat$Lat), to =max(myDat$Lat), by =1) # get a range of latitudes for making predictionsforSex<-unique(myDat$Sex) # get every level of my Sex covariateforVis<-expand.grid(Lat=forLat, Sex=forSex) # create a data frame with all combinations of covariates# Get your model fit estimates at each value of your covariatesmodFit<-predict(bestMod, # the modelnewdata = forVis, # the covariate valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frameggplot()+geom_point(data = myDat, # datamapping =aes(x = Lat, y = delMass, col = Sex))+# add data to your plotgeom_ribbon(data = forVis, mapping =aes(x = Lat, ymin = Lower, ymax = Upper, fill = Sex), alpha =0.5)+# add the uncertainty to your plotgeom_line(data = forVis, mapping =aes(x = Lat, y = Fit, col = Sex))+# add the model fit to your plotylab("Change in body mass (kg per year)")+# change y-axis labelxlab("Latitude (˚N)")+# change x-axis labellabs(caption ="Figure 4: Modelled effects (lines, mean +/- standard error) of latitude and sex (colours) on growth rate \nalong with observations (points)")+theme_bw()+# change ggplot themetheme(plot.caption =element_text(hjust=0)) # move figure legend (caption) to left alignment. Use hjust = 0.5 to align in the center.
# * Report how well your model explains your response (explained deviance)r.squaredLR(bestMod) # estimates the likelihood ratio R^2. Could also estimate a traditional R^2 with 1-summary(bestMod)$deviance/summary(bestMod)$null.deviance here as the error distribution assumption is normal and the shape assumption is linear, but the likelihood ratio R^2 function is generally applicable to many error distribution assumptions and equivalent to the traditional R^2 when the error distribution assumption is normal and the shape assumption is linear.
# * report how important each model term (covariate fixed effects and any interactions) is in explaining the deviation in your response. sw(dredgeOut)
Sex Lat Lat:Sex
Sum of weights: 1.00 1.00 0.76
N containing models: 3 3 1
The write-up:
Use the results of the code above to report your hypothesis testing results. Include:
a description of the terms (covariate fixed effects and any interactions) that are in your best-specified model
compare these terms to your starting model - are all terms in your starting model in you best-specified model?
report your modelled effects by reporting your coefficients. As we will discuss this week in class, report:
your coefficient estimates with uncertainty
evidence that these coefficient estimates are different than zero
evidence that these coefficient estimates are different than each other across the levels of your categorical covariate (if you have one)
report your modelled effects by plotting. Remember to include your model predictions, uncertainty as well as your observations. Include units on your axes and a figure number and legend.
report your estimate of the % of deviance in your response that your model explains
consider what might be explaining the remaining (unexplained) deviation in your response
report how important each model term (covariate fixed effects and any interactions) is in explaining the deviation in your response.
For example,
My best-specified model tells me that the growth varies with latitude and sex, and that the effect of latitude on growth varies with sex. The terms in my best-specified model are the same as those in my starting model indicating that there is evidence that the main effects of latitude and sex as well as the interaction are explaining variability in moose growth.
Growth rate is higher for males than females: the predicted growth when latitude is 61.6˚N (the mean of the latitudinal range) is 41.8 ± 0.6 kg year-1 for females and 53.8 ± 0.6 kg year-1 for males.40 Both these coefficients are significantly different than zero (t-test, P < 0.0001 for males and females)^[note that we don’t need to add an adjustment for multiple comparisons here as there are only two levels in our categorical covariate).
The coefficient (slope) for the effect of a change of latitude on the growth rate is 0.70 ± 0.18 kg year-1 ˚N-1 for females and 1.29 ± 0.18 kg year-1 ˚N-1 for males. These slopes are significantly different from zero (t-test, P = 0.001 for females and P < 0.0001 for males) and significantly different from one another (t-test, P = 0.035). These coefficients show that growth rates increase with latitude for both sexes, but that the effect of latitude on growth is higher for male vs. female moose.
Figure 4 shows the modelled effects of latitude and sex on growth rate along with my observations.
Together the effects of latitude and sex explain 93% of the deviance in growth rate (Likelihood ratio R2).41
The remaining unexplained deviance may be due to other factors affecting growth rate such as winter harshness, food availability, etc. Note that I was not able to include minimum winter temperature in my hypothesis test due to high collinearity between latitude and minimum winter temperature. Therefore, I do not know if the measured latitudinal effect on growth might mechanistically be due to winter harshness (here estimated as minimum winter temperature). This could be tested if I was able to expand my data set to include sites where latitude and minimum winter temperature were less correlated.
Based on the sum of Akaike weights, latitude and sex are equally important in explaining the deviance in growth rate as they appear in all models with Akaike weights > 0. The interaction between latitude and sex is slightly less important appearing only in the best-specified model with an Akaike weight of 0.76.
6.9 Predicting
The task:
Here you may use your model to make predictions (remembering prediction limits). For the assignments and exam, I’ll explicitly ask you to make a prediction if we want to see one.
If you want to make a prediction (e.g. based on your best-specified model, what is your predicted mean annual growth of a female moose living in Sicily, Italy?), in this section you want to report your prediction results and give any limitations to the prediction.
The code:
Code supporting your tasks in this section will help you communicate what prediction you made and the results, e.g.:
# Based on your best-specified model, what is your predicted mean annual growth of an adult female moose living in Sicily, Italy?## Looking on the internet, we see that Sicily, IT is at 37.6˚N.predict(bestMod, # our modelnewdata =data.frame(Sex ="Females", Lat =37.6), # values of the covariates at which to make the predictionse.fit=TRUE, # include an estimate of error around the predictiontype="response") # make sure the prediction is on the response scale
Use the code in the previous section to present the results of your prediction. Include
the prediction and an estimate of uncertainty
any perceived prediction limits.
For example,
I used my best-specified model to estimate the mean annual growth rate of adult female moose living in Sicily, Italy (~ 37.6˚N) as 24.9 +/- 4.4 kg per year (mean +/- standard error). While this is an estimate consistent with my model, it likely is unrealistic as the distribution of moose does not currently extend as far south as Sicily. It is likely there are limitations to their dispersal to or survival in this area and my model may not be valid for latitudes so far south.
- Reporting your best-specified model - Reporting your modelled patterns - Plotting model fit, uncertainty and observations [[Visualizing]] - Reporting effect sizes - Report the coefficients - [[what are model coefficients]] - What are the coefficient estimates (with uncertainty)? - Is each coefficient significantly different from zero? - Are the coefficients significantly different than one another? - Reporting how well your model explains your response (explained variation) - Report the relative importance of each term in your model
cross validation
visreg vs. ggeffects
three ways to report what your model is saying about your hypothesis
- visualize
- report coefficients
- give an example (response predicted under conditions)
- model coefficients are model effects
- the effect on Resp with a change in Cov
- Continuous
- coefficient is the change in Resp with a unit change in Cov
- what does it mean if coef is positive
- what does it mean if coef is negative
- what is error on the coef?
- what does it mean when a coef is no different than 0? - why that is the science
- how to test if two coefficients are similar
- Categorical
- coefficient is the change in Resp with a change from one level of Cat to another
- what does it mean if coef is positive
- what does it mean if coef is negative
- what is error on the coef?
- what does it mean when a coef is no different than 0? - why that is the science
- how to test if two coefficients are similar
- Show linear model equation and how it looks for Cat, Cont with coefficients
- Describe generically and then define for each link function
report effect size and precision, not p-value
Emphasise effect sizes to replace statistical significance with ecological relevance
research question should follow the FINER criteria
refinements of your hypothesis after you collect the data will result in a hypothesis that reflects the sample, not the population
exploratory or hypothesis generating analysis
use representative (unbiased) sampling
to report results of hypothesis testing, report effect sizes, confidence intervals, interpretation of effect and confidence intervals, and can include p-values
p-values are the probability that we would obtain a result equal to or more extreme than that observed given a null hypothesis that is true
p-values tell us more about how compatible the data are with the model vs. evidence of underlying effects
7 A statistical modelling example
It will help our discussions to have an example of a best-specified model. For our example, let’s assume that your hypothesis is that your response variable Resp is explained by a categorical covariate (Cat), a continuous covariate (Cont) as well as the interaction between the two. You can communicate this as:
Resp ~ Cat + Cont + Cat:Cont
Let’s assume you tested this hypothesis by fitting a model with a Gamma error distribution assumption (to reflect the nature of your response variable) and linear shape assumption (to reflect the relationship among each covariate and your response) using a GLM:
You considered possible covariate collinearity42 and observation independence43 and how it might be affecting your starting model. You then assessed your starting model to make sure the error distribution and shape assumptions were appropriate using the DHARMa method:
library(DHARMa)simulationOutput <-simulateResiduals(fittedModel = startMod, n =250) # simulate data from our model n timesplot(simulationOutput, asFactor=FALSE) # compare simulated data to our observations
plot(simulationOutput, # compare simulated data to form=myDat$Cont, # our observationsasFactor=FALSE) # whether the variable plotted is a factor
plot(simulationOutput, # compare simulated data to form=myDat$Cat, # our observationsasFactor=TRUE) # whether the variable plotted is a factor
All looks fine - the distribution of the expected and observed residuals are similar, and there is no pattern in the residuals when plotted with fitted values or any of the covariates.
You then used your starting model to test your hypothesis by model selection using the dredge() function in the MuMIn package:
library(MuMIn)options(na.action ="na.fail") # needed for dredge() function to prevent illegal model comparisons(dredgeOut<-dredge(startMod, extra ="R^2")) # fit and compare a model set representing all possible predictor combinations
You find two models have something to say about your hypothesis but the top model (model #8) is greater than 2 AICc lower than the second model (model #4). You choose your best-specified model as:
bestMod<-(eval(attr(dredgeOut, "model.calls")$`8`)) # extract model #8 from dredge tablebestMod
Here we’ll use this example best-specified model to visualize and report what it is saying about your hypothesis. We’ll cover other examples in class as well.
8 Visualizing
Visualizing your model means visualizing what your model is saying about your research hypothesis. We have practiced visualizing your fitted statistical model but now let’s focus on visualizing to communicate the effects described in your model.
Here, you want to visualize how each covariate is affecting the response by drawing how fitted values (the estimates of your response made by the model) change with your covariates, along with a measure of uncertainty around those fitted values. You also want to include your data (actual observations of your response and covariates) so you can start to communicate how well your model captures your hypothesis, and what amount of deviance in your response is explained by your model.
There are a lot of different R packages that make it easy to quickly visualize your model. We’ll focus on two methods that will allow you to make quick visualizations of your model and/or customize the figures as you would like.
8.1 Visualizing with the visreg package
The first method uses the visreg package and the visreg() function which we’ve practiced before. Using your example bestMod, here are a few visualizations using the visreg package:
Visualizing with a continuous covariate on the x-axis: Notice how you can include the gg = TRUE argument to plot this as a ggplot type figure. This allows you to add your data onto the visualization of your model.
library(visreg) # load visreg packagelibrary(ggplot2) # load ggplot2visreg(bestMod, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cont", # covariate on x-axisby ="Cat", # covariate plotted as colour#breaks = 3, # if you want to control how the colour covariate is plotted#cond = , # if you want to include a 4th covariateoverlay =TRUE, # to plot as overlay or panels rug =FALSE, # to include a ruggg =TRUE)+# to plot as a ggplotgeom_point(data = myDat, # datamapping =aes(x = Cont, y = Resp, col = Cat))+# add data to your plotylim(0, 60)+# adjust the y-axis unitsylab("Response, (units)")+# change y-axis labelxlab("Cont, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
Visualizing with a categorical covariate on the x-axis: You can also plot the same model with the categorical covariate on the x-axis. Notice how the continuous covariate is represented by different levels in the colours on the plot. Here you’ve asked for three levels with breaks = 3.
visreg(bestMod, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cat", # covariate on x-axisby ="Cont", # covariate plotted as colourbreaks =3, # if you want to control how the colour covariate is plotted#cond = , # if you want to include a 4th covariateoverlay =TRUE, # to plot as overlay or panels rug =FALSE, # to include a ruggg =TRUE)+# to plot as a ggplotylab("Response, (units)")+# change y-axis labelxlab("Cat, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
You can also specify at which levels the breaks should occur with the breaks = ... argument. For example, you can ask visreg to plot the modelled effects when Cont = 40 and Cont = 50:
visreg(bestMod, # model to visualizescale ="response", # plot on the scale of the responsexvar ="Cat", # covariate on x-axisby ="Cont", # covariate plotted as colourbreaks =c(40,50), # if you want to control how the colour covariate is plotted#cond = , # if you want to include a 4th covariateoverlay =TRUE, # to plot as overlay or panels rug =FALSE, # to include a ruggg =TRUE)+# to plot as a ggplotylab("Response, (units)")+# change y-axis labelxlab("Cat, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
Unfortunately due to limitations of the visreg package, you can not easily add our data onto plots where the x-axis is a categorical covariate. But that’s ok, because there are other options…
8.2 Visualizing “by hand”:
We will also practice how to visualize your model “by hand”. Here, “by hand” is a bit of a silly description as R will be doing the work for you. What I mean by “by hand” is that you will be building the plot yourselves by querying your model object.
Visualizing with a continuous covariate on the x-axis:. To plot by hand, you will first make a data frame containing the value of your covariates at which you want to plot effects on the response:
# Set up your covariates for the visualized fitforCont<-seq(from =min(myDat$Cont), to =max(myDat$Cont), by =0.1) # e.g. a sequence of Cont valuesforCat<-unique(myDat$Cat) # every value of your categorical covariate# create a data frame with your covariatesforVis<-expand.grid(Cont=forCont, Cat=forCat) # expand.grid() function makes sure you have all combinations of covariates
Next, you will use the predict() function44 to get the model estimates of your response variable at those values of your covariates:
# Get your model fit estimates at each value of your covariatesmodFit<-predict(bestMod, # the modelnewdata = forVis, # the covariate valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame
Finally, you will use these model estimates to make your plot:
ggplot()+geom_point(data = myDat, # datamapping =aes(x = Cont, y = Resp, col = Cat))+# add data to your plotgeom_ribbon(data = forVis, mapping =aes(x = Cont, ymin = Lower, ymax = Upper, fill = Cat), alpha =0.5)+# add the uncertainty to your plotgeom_line(data = forVis, mapping =aes(x = Cont, y = Fit, col = Cat))+# add the model fit to your plotylim(0, 60)+# adjust the y-axis unitsylab("Response, (units)")+# change y-axis labelxlab("Cont, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
Visualizing with a categorical covariate on the x-axis:. The procedure for visualizing when you have a categorical covariate on the x-axis is similar: 1) choose the values of your covariates at which to make a prediction, 2) use predict() to use your model to estimate your response variable at those values of your covariate, and 3) use the model estimates to plot your model fit:
# Set up your covariates for the visualized fitforCont<-c(40, 50 ,60) # e.g. Cont valuesforCat<-unique(myDat$Cat) # every value of your categorical covariate# create a data frame with your covariatesforVis<-expand.grid(Cont=forCont, Cat=forCat) # expand.grid() function makes sure you have all combinations of covariates# Get your model fit estimates at each value of your covariatesmodFit<-predict(bestMod, # the modelnewdata = forVis, # the covariate valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frameggplot()+geom_jitter(data = myDat, # datamapping =aes(x = Cat, y = Resp))+# add data to your plotgeom_errorbar(data = forVis, mapping =aes(x = Cat, y = Fit, col = Cont, ymin = Lower, ymax = Upper))+# add the uncertainty to your plotylim(0, 60)+# adjust the y-axis unitsylab("Response, (units)")+# change y-axis labelxlab("Cat, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
8.3 Plotting in 3D
You might notice the plots above are communicating 3-dimensions (one response + two covariates) in a 2-dimensional plot. There are other ways of making 3-dimensional plots in R, e.g. with the visreg package using the visreg2d() function in the visreg package:
visreg2d(fit = bestMod, # your modelxvar ="Cont", # what to plot on the x-axisyvar ="Cat", # what to plot on the y-axisscale ="response") # make sure fitted values (colours) are on the scale of the response
or “by hand” using the geom_raster() function in the ggplot2 package:
# Set up your covariates for the visualized fitforCont<-seq(from =min(myDat$Cont), to =max(myDat$Cont), by =0.1) # e.g. a sequence of Cont valuesforCat<-unique(myDat$Cat) # every value of your categorical covariate# create a data frame with your covariatesforVis<-expand.grid(Cont=forCont, Cat=forCat) # expand.grid() function makes sure you have all combinations of covariates# Get your model fit estimates at each value of your covariatesmodFit<-predict(bestMod, # the modelnewdata = forVis, # the covariate valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data frame# create your plotggplot() +# start your ggplotgeom_raster(data = forVis, aes(x = Cont, y = Cat, fill = Fit))+# add the 3 dimensions as a rastergeom_point(data = myDat, aes(x = Cont, y = Cat, colour = Resp)) # add your data
EXTRA: As you can see, these plots represent the 3rd dimension by using colour. We can also make actual 3 dimensional plots in R with the plotly package45. These plots are interactive which makes them more useful than static 3d plots. Use your mouse to move the plot around!
library(plotly) # load the plotly packageplot_ly(data = forVis, # the data with your model predictions (made above)x =~Cont, # what to plot on the x axisy =~Cat, # what to plot on the y axisz =~Fit, # what to plot on the z axistype="scatter3d", # type of plotmode="markers") %>%# type of plotadd_markers(data = myDat, x =~Cont, y =~Cat, z =~Resp) # add your data
9 Reporting
To report your model, you want to communicate a number of things:
* Reporting your best-specified model(s). * Reporting the modelled effects of your covariates on your response.
* Reporting how well your model explains your response. * Reporting the importance of each covariate in explaining your response.
9.1 Reporting your best-specified model(s)
Reporting your best-specified model means reporting the terms - the covariates and any interactions - that are in your model. It’s good practice to present the model along with the results from the model selection. In this way, you can include multiple best-specified models if there is evidence that more than one might be useful.
For example, with your results from model selection being:
Depending on your hypothesis and results, you may want to present all models within ∆AICc < 2 of the best-specified model, or all models with any Akaike weight (as above), or simply all models. We’ll discuss ways of making these choices in class.
9.2 Reporting the modelled effects of your covariates on your response.
You want to report how the covariates are affecting your response. e.g. does a change in the covariate lead to an increase in the response? or a decrease? or is the effect non-linear (for continuous covariates)?
The easiest way to do this is by visualization, as you saw above, but you also want to report the coefficients with their uncertainty:
Note that for categorical covariates with more than two levels (e.g. an example with a covariate consisting of Species A, B and C), you’ll want to take an extra step to determine if the effect of each level of the covariate is similar (e.g. is the effect of Species A on the response the same as the effect on Species B?). We’ll cover how to do this in Week 12.
9.3 Reporting how well your model explains your response
If you will recall, your whole motivation for pursuing statistical modelling was to explain variation in your response. Thus, it is important that you quantify how much variation in your response you are able to explain by your model.
Note that we will discuss this as “explained deviance” rather than “explained variation”. This is because the term “variance” is associated with models where the error distribution assumption is normal, whereas deviance is a more universal term.
When you have a normal error distribution assumption and linear shape assumption, you can capture the amount of explained deviance simply by comparing the variance in your response (null variation) vs. the variance in your model residuals (residual variation) as the \(R^2\):
From this equation, you can see how, if your model is able to explain all the variation in the response, the residual variation will be zero, and \(R^2 = 1\). Alternatively, if the model explains no variation in the response the residual variation equals the null variation and \(R^2 = 0\).
For models with other error distribution and shape assumptions, you need another way of estimating the goodness of fit of your model. You can do this through a pseudo \(R^2\).
One useful pseudo \(R^2\) is called the Likelihood Ratio \(R^2\) or \(R^2_{LR}\). The \(R^2_{LR}\) compares the likelihood of your best-specified model to the likelihood of the null model:
where \(n\) is the number of observations, \(log𝓛(model)\) is the log-likelihood of your model, and \(log𝓛(null)\) is the log-likelihood of the null model. The \(R^2_{LR}\) is the type of pseudo \(R^2\) that shows up in your dredge() output when you add extra = "R^2" to the dredge() call. You can calculate \(R^2_{LR}\) by hand, read it from our dredge() output, or estimate it using r.squaredLR() from the MuMIn package:
Note here that two values of \(R^2_{LR}\) are reported. The adjusted pseudo \(R^2\) given here under attr(,"adj.r.squared") has been scaled so that \(R^2_{LR}\) can reach a maximum of 1, to be equivalent to a regular \(R^2\).
One nice feature of the \(R^2_{LR}\) is that it is equivalent to the regular \(R^2\) when our model assumes a normal error distribution assumption and linear shape assumption, so you can use\(R^2_{LR}\) for any of the models we’re discussing in class.
check RR2 package
9.4 Reporting the importance of each covariate in explaining your response.
Finally, you want to report how relatively important each term is in your model in regards to explaining deviance in your response. This is also called “partitioning the explained deviance” to the covariates.
To get an estimate of how much deviance in your response one particular covariate explains, you may be tempted to compate the explained deviance (\(R^2\)) estimates of models fit to data with and without that particular covariate. Let’s try an example using your example statistical model:
If you want to get an estimate as to how much response deviance the Cont covariate explains, you might compare the pseudo \(R^2\) of a model with and without the Cont covariate, e.g. comparing model #4 (that includes Cont) and model #2 (that doesn’t include Cont):
# deviance explained by Cont:R2.mod4 <- (dredgeOut$`R^2`[2]) # model #4 is the second row in dredgeOutR2.mod2 <- (dredgeOut$`R^2`[3]) # model #2 is the third row in dredgeOutR2.mod4 - R2.mod2 # find estimated contribution of Cont to explained deviance
[1] 0.1731384
But what if you chose to compare model #3 (that includes Cont) and model #1 (that doesn’t include Cont):
# deviance explained by Cont:R2.mod3 <- (dredgeOut$`R^2`[4]) # model #3 is the fourth row in dredgeOutR2.mod1 <- (dredgeOut$`R^2`[5]) # model #1 is the fifth row in dredgeOutR2.mod3 - R2.mod1 # find estimated contribution of Cont to explained deviance
[1] 0.372649
Quite a different answer! Your estimates of the contribution of Cont to explaining the response deviation don’t agree because of collinearity among our covariates46. We’ll discuss this more in class.
One option for partitioning the explained deviance when you have collinearity among your covariates is hierarchical partitioning. Hierarchical partitioning estimates the average independent contribution of each covariate to the total explained variance by considering all models in the candidate model set47. This method is beyond the scope of the course but see the rdacca.hp package for an example of how to do this.
Another method for estimating the importance of each term (covariate or interaction) in your model is by again looking at the candidate model set ranking made by dredge(). Here you can measure the importance of a covariate by summing up the Akaike weights for any model that includes a particular covariate. The Akaike weight of a model compares the likelihood of the model scaled to the total likelihood of all the models in the candidate model set. The sum of Akaike weights for models including each covariate tells us how important the covariate is in explaining the deviance in your response. You can calculate the sum of Akaike weights with the sw() function in the MuMIn package:
sw(dredgeOut)
Cat Cont Cat:Cont
Sum of weights: 1.00 1.00 0.82
N containing models: 3 3 1
Here we can see that Cat and Cont are equally important in explaining the deviance in Resp (they appear in all models that have any Akaike weight), while the interaction term between Cat and Cont is less important (only appearing in one model with Akaike weight, though this is the top model).
Check relaimpo
10 Predicting
Now that you have your best-specified model, you may want to use it to make a prediction, e.g. “what would I expect the value of my response to be at a different value of my covariate?” This may include estimating your response at a different place, time or environmental (covariate) condition.
10.1 First, consider the limits to making predictions
The first step when considering using your model to make predictions is to decide if this is appropriate. By using your model to make a prediction, you are assuming that the same processes that govern the relationship between your covariate(s) and response will hold in the new place, time or environmental condition. For example, you might have used a linear shape assumption that proved appropriate for modelling how growth varies due to temperature, but now you want to estimate growth at an even higher temperature than you tested. Is this appropriate? Would it be possible that growth starts to respond non-linearly with temperature?
10.2 How to use your model to make predictions
Once you are convinced that making a prediction with your model is useful, you can use the predict() function to predict the value of your response by handing it values of each covariate at which you would like a response prediction. For example:
predict(object = bestMod, # your modelnewdata =data.frame(Cont =57, Cat ="Treatment"), # the values of the covariates at which to make the predictiontype ="response", # to make the prediction on the response scale.se.fit =TRUE) # to include a measure of uncertainty around the prediction
So with Cont = 57 and Cat = Treatment, we would expect Resp to be 34.6 ± 1.1. This can be shown on our visualization here (prediction in black):
11 Reporting models with categorical covariates
In this week’s notes, we’ll work through three example best-specified models to report our model results, and in some cases make predictions.
Note that each example is testing a different hypothesis by fitting a different model to a different data-set. The three examples are not related to one another.
Note that I’m not going through all the steps that got us to these example best-specified models, but assume we went through the previous steps (Response, Covariates, Hypothesis, Starting model, Assessing Fit and Hypothesis Test) as we practiced in class to choose each best-specified model.
11.1 Example 1: Resp1 ~ Cat1
For example #1, assume your best-specified model shows that your response variable (Resp1) is explained by a categorical covariate (Cat1):
Resp1 ~ Cat1
Your best-specified model for example #1 is in an object called bestMod1:
bestMod1
Call: glm(formula = Resp1 ~ Cat1 + 1, family = gaussian(link = "identity"),
data = myDat1)
Coefficients:
(Intercept) Cat1Sp2 Cat1Sp3
22.500 4.194 -5.889
Degrees of Freedom: 99 Total (i.e. Null); 97 Residual
Null Deviance: 8882
Residual Deviance: 7037 AIC: 717.2
that was fit to data in myDat1:
summary(myDat1)
Cat1 Resp1
Sp1:28 Min. : 1.00
Sp2:36 1st Qu.:15.00
Sp3:36 Median :22.00
Mean :21.89
3rd Qu.:30.00
Max. :40.00
and the dredge() table you used to pick your bestMod1 is in dredgeOut1
dredgeOut1
Global model call: glm(formula = Resp1 ~ Cat1, family = gaussian(link = "identity"),
data = myDat1)
---
Model selection table
(Intrc) Cat1 R^2 df logLik AICc delta weight
2 22.50 + 0.2077 4 -354.584 717.6 0.00 1
1 21.89 0.0000 2 -366.223 736.6 18.98 0
Models ranked by AICc(x)
11.1.1 Reporting your statistical modelling results - Example 1
Recall that reporting your statistical modelling results means:
Reporting your best-specified model(s)
Reporting your modelled effects
Reporting how well your model explains your response
Reporting how important each of your covariates is in explaining your response
Reporting best-specified model(s) - Example 1
Reporting your best-specified model means reporting the terms - the covariates and any interactions - that are in your model.
For Example 1, you will report that your best-specified model is Resp1 ~ Cat1, i.e. that there is evidence that Cat1 explains variability in Resp1. Remember from last week that you can also use the output from dredgeOut1 to report evidence for how you picked the best-specified model, e.g. if you want to report more than one model.
Global model call: glm(formula = Resp1 ~ Cat1, family = gaussian(link = "identity"),
data = myDat1)
---
Model selection table
(Intrc) Cat1 R^2 df logLik AICc delta weight
2 22.50 + 0.2077 4 -354.584 717.6 0.00 1
1 21.89 0.0000 2 -366.223 736.6 18.98 0
Models ranked by AICc(x)
Reporting your modelled effects - Example 1
Reporting the coefficients
If you remember from last week, you reported your modelled effects by reporting your model coefficients as your coefficients tell you about the effect your covariate has on your response.
With a continuous covariate, the coefficient is the slope of the line (on the linked scale) showing the amount of change in the response that is caused by a change in your continuous covariate.
With a categorical covariate, the coefficient is an intercept, or the model prediction of the response at that level of the categorical covariate (e.g. the modelled Resp1 when Cat1 is Sp1, Sp2 and Sp3). For categorical covariates, there is a coefficient estimated for each level of the covariate that can be thought of as an intercept. So as Cat1 in Example 1 has three levels (Sp1, Sp2 and Sp3), there will be three coefficients estimated (one for each).
Here we will focus on three things:
What are your coefficient estimates (with uncertainty)?
Is each coefficient estimate different from zero?
Are the coefficient estimates different than one another across categorical levels?
What are your coefficient estimates (with uncertainty)?
Last week, you found your coefficients in the “Estimate” column of the summary() output of your model:
coef(summary(bestMod1)) # extract the coefficients from summary()
Interpreting the coefficients from this output takes practice - especially for categorical covariates because of the way that R treats categorical covariates in regression. With R’s “dummy coding”, one level of the covariate (here Sp1) is incorporated into the intercept estimate 22.5 as the reference level. The other coefficients in the Estimate column represent the change in modelled response when you move from the reference level (Sp1) to another level.
The modelled response when Cat1 = Sp1 is 22.5 units.
The modelled response when Cat1 = Sp2 is 4.2 units higher than this reference level = 22.5 + 4.2 = 26.7 units.
The modelled Resp1 when Cat1 = Sp3 is -5.9 units lower than the reference level = 22.5 + -5.9 = 16.6 units.48
So all the information we need is in this summary() output, but not easy to see immediately. An easier way is to use the emmeans package which helps us by reporting the coefficients directly. “emmeans” is the estimated marginal means. The estimated marginal mean is the model prediction for the response when the covariates are held at particular values. If the value of the covariate is not defined, the response prediction is the average response prediction over all values in that covariate.
In our case, we use the emmeans package to get the mean value of the response at each level of the covariate. For categorical covariates, we do this with the emmeans() function:
library(emmeans) # load the emmeans packageemmeans(object = bestMod1, # your modelspecs =~ Cat1, # your covariatestype ="response") # report coefficients on the response scale
So now we have a modelled value of our response for each level of our categorical covariate.
When Cat1 is Sp1, Resp1 is estimated to be 22.5 ± 1.6 units.
When Cat1 is Sp2, Resp1 is estimated to be 26.7 ± 1.4 units.
When Cat1 is Sp3, Resp1 is estimated to be 16.6 ± 1.4 units.
Note that this is the same information in the summary() output just easier to read.49
Note also that two types of uncertainty are measured here. SE stands for the standard error around the prediction, and is a measure of uncertainty of the average modelled effect. The lower.CL and upper.CL represent the 95% confidence limits of the prediction - so if I observed a new Resp1 at a particular Cat1, there would be a 95% chance it would lie between the bounds of the confidence limits.
Is each coefficient estimate different from zero?
Cat1 is in your best-specified model, so you know that there is evidence it has an effect50 on your response (Resp1). But is there an effect on Resp1 for all levels of Cat1 - are any of the coefficients statistically similar to 0?
forComp <-emmeans(object = bestMod1, # your modelspecs = pairwise ~ Cat1, # your covariates, with the request for a pairwise testtype ="response") # report coefficients on the response scaletest(forComp)$emmeans # shows test of coefficients = 0
Here we see that the coefficients for all levels are significantly different than zero (t-test, P < 0.0001.)
Are the coefficient estimates different than one another across categorical levels?
Finally, we ask if all levels of Cat1 affect Resp1 in the same way - i.e. are the coefficients across factor levels significantly different from one another?
To get evidence about how each level in your categorical covariate affects your response, you need to test which effects differ from one another using multiple comparisons51, i.e. you will compare the modelled effect of each level of your categorical covariate vs. each other level of your categorical covariate to determine which are different from each other (called pairwise testing). You will do this by testing the null hypothesis that the modelled effects are similar to one another, typically rejecting the null hypothesis when p < 0.05. Remember that the p-value is the probability that we observe a value at least as big as the one we observed even though our null hypothesis is true. In this case, we are looking at the difference between coefficients estimated for two levels of our covariate. The p-value is the probability of getting a difference at least as big as the one we observed even though there is actually no difference between the coefficients (the null hypothesis is true).
A couple of things to note about multiple comparison testing:
Multiple comparison testing on a categorical covariate should only be done after your hypothesis testing has given you evidence that the covariate has an effect on your response. That is, you should only do a multiple comparison test on a covariate if that covariate is in your best-specified model. As this is a test done after your hypothesis testing, it is called a post-hoc test.
Multiple comparison testing can be a problem because you are essentially repeating a hypothesis test many times on the same data (i.e. are the effects of Sp1 different than those of Sp2? are the effects of Sp1 different than those of Sp3? are the effects of Sp2 different than those of Sp3?…). These repeated tests mean there is a high chance that you will find a difference in one test purely due to random chance, vs. due to there being an actual difference. To account for this, the multiple comparison tests you will perform have been formulated to correct for this increased error.
Multiple comparison testing is very simple with the emmeans package. It just requires us to add that we would like to have multiple comparisons testing in our call by adding pairwise to the specs = argument:
forComp <-emmeans(object = bestMod1, # your modelspecs = pairwise ~ Cat1, # your covariates, with the request for a pairwise testtype ="response") # report coefficients on the response scaletest(forComp)
$emmeans
Cat1 emmean SE df t.ratio p.value
Sp1 22.5 1.61 97 13.978 <.0001
Sp2 26.7 1.42 97 18.804 <.0001
Sp3 16.6 1.42 97 11.701 <.0001
$contrasts
contrast estimate SE df t.ratio p.value
Sp1 - Sp2 -4.19 2.15 97 -1.954 0.1292
Sp1 - Sp3 5.89 2.15 97 2.744 0.0196
Sp2 - Sp3 10.08 2.01 97 5.023 <.0001
P value adjustment: tukey method for comparing a family of 3 estimates
You can see that you get two tables from this new emmeans() call. The first table ($emmeans) is the same data you saw before - the modelled coefficients (intercepts) for each level in Cat1. The second table ($contrasts) shows the results of the multiple comparison (pairwise) testing. The values in the p.value column tell us the results of the hypothesis test comparing the coefficients between the two levels. For example, for Sp1 vs. Sp2, there is a 13% (p = 0.13) probability of getting a difference in coefficients at least as big as 22.5 - 26.7 = -4.2 even though the null hypothesis (no difference) is true. This value 13% (p = 0.13) is too big (i.e. bigger than our threshold of 5% or p = 0.05) for us to believe we have evidence the coefficients are different from one another.
Based on a threshold p-value of 0.05, we can see that:
The value of Resp1 when Cat1 is Sp1 is not statistically different than that when Cat1 is Sp2 as p = 0.13 is greater than p = 0.05.
The value of Resp1 when Cat1 is Sp1 is different (higher) than that when Cat1 is Sp3 as p = 0.02 is smaller than p = 0.05. The value of Resp1 when Cat1 is Sp2 is different (higher) than that when Cat1 is Sp3 as p < 0.000152 is smaller than p = 0.05.53.
Note that the p-values have been adjusted via the Tukey method which adjusts the difference that the two coefficients need to have to allow for the fact that we’re making multiple comparisons.54
Note that you can also get the results from emmeans visually with
plot(forComp,comparisons =TRUE)
Plotting modelled effects
And this brings us to another important aspect of reporting: visualizing your modelled effects. As mentioned last class, you want to include i) your model fit, ii) uncertainty around that fit, and iii) your observations on your plot. You can do this “by hand” with:
# Set up your covariates for the visualized fitforCat1<-unique(myDat1$Cat1) # every value of your categorical covariate# create a data frame with your covariatesforVis<-expand.grid(Cat1=forCat1) # expand.grid() function makes sure you have all combinations of covariates# Get your model fit estimates at each value of your covariatesmodFit<-predict(bestMod1, # the modelnewdata = forVis, # the covariate valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data framelibrary(ggplot2)ggplot() +# start ggplotgeom_point(data = myDat1, # add observations to your plotmapping =aes(x = Cat1, y = Resp1), position=position_jitter(width=0.1)) +# control position of data points so they are easier to see on the plotgeom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat1, y = Fit, ymin = Lower, ymax = Upper),size=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat1, y = Fit), shape =22, # shape of pointsize =3, # size of pointfill ="white", # fill colour for plotcol ='black') +# outline colour for plotylab("Resp1, (units)")+# change y-axis labelxlab("Cat1, (units)")+# change x-axis labeltheme_bw() # change ggplot theme
Reporting how well your model explains your response - Example 1
As a reminder from last week, you can estimate the deviance explained in your response by using a pseudo \(R^2\) called the Likelihood Ratio \(R^2\) or \(R^2_{LR}\).55. You can calculate \(R^2_{LR}\) by hand, read it from our dredge() output, or estimate it using r.squaredLR() from the MuMIn package:
Note here that two values of \(R^2_{LR}\) are reported. The adjusted pseudo \(R^2\) given here under attr(,"adj.r.squared") has been scaled so that \(R^2_{LR}\) can reach a maximum of 1, to be equivalent to a regular \(R^2\).
So you estimate that ~ 21% of the deviance in Resp1 is explained by Cat1 via bestMod1.
One nice feature of the \(R^2_{LR}\) is that it is equivalent to the regular \(R^2\) when your model assumes a normal error distribution assumption and linear shape assumption, so you can use\(R^2_{LR}\) for any of the models we’re discussing in class. Let’s check this here by comparing our \(R^2_{LR}\) to that from true \(R^2\) we calculated in Week 9:
\(R^2_{LR}\) and \(R^2\) are identical because you have a normal error distribution assumption and linear shape assumption.
Reporting covariate importance - Example 1
With Example 1, you have only one covariate (Cat1) and so this covariate is responsible for explaining all of the variability in your response (Resp1). You can see that it appears in all models with any weight with your sw() function from the MuMIn package.
dredgeOut1
Global model call: glm(formula = Resp1 ~ Cat1, family = gaussian(link = "identity"),
data = myDat1)
---
Model selection table
(Intrc) Cat1 R^2 df logLik AICc delta weight
2 22.50 + 0.2077 4 -354.584 717.6 0.00 1
1 21.89 0.0000 2 -366.223 736.6 18.98 0
Models ranked by AICc(x)
sw(dredgeOut1)
Cat1
Sum of weights: 1
N containing models: 1
11.1.2 Predicting with Example 1
With Example 1, we have only a categorical covariate. In such cases, it doesn’t make sense to use your model to make a prediction as you can’t make a prediction for a level not already in your covariate, e.g. there’s no way to predict Resp1 if Cat1 was Sp5 or Sp89.
So if you have only categorical covariates, you aren’t able to use your model for prediction. Coming up (Example 3) we’ll look at an example with a combination of continuous and categorical covariates and, in these cases, prediction is possible.
11.2 Example 2: Resp2 ~ Cat2 + Cat3 + Cat2:Cat3
For Example #2, assume your best-specified model is that your response variable (Resp2) is explained by two categorical covariates (Cat2 & Cat3) as well as the interaction between the covariates:
Resp2 ~ Cat2 + Cat3 + Cat2:Cat3
Your best-specified model for example #2 is in an object called bestMod2:
For Example 2, you will report that your best-specified model is Resp2 ~ Cat2 + Cat3 + Cat2:Cat3, i.e. that there is evidence that Cat2 and Cat3 explain variability in Resp2, and that there is an interaction between Cat2 and Cat3 - i.e. the effect of Cat2 on Resp2 depends on the value of Cat3. Remember that you can also use the output from dredgeOut2 to report evidence for how you picked the best-specified model, e.g. if you want to report more than one model.
What are your coefficient estimates (with uncertainty)?
Again, for categorical covariates, there is a coefficient estimated for each level of the covariate. If there is one or more interactions among covariates, there will be one coefficient for each combination of levels among covariates. Let’s look at the summary() output of your model to understand better:
coef(summary(bestMod2)) # extract the coefficients from summary()
Comparing this output to that of Example 1 shows many more estimated coefficients for Example 2. This is because we have one coefficient for each level of each covariate, as well as coefficients for each combination of levels of the covariates.
Again, it takes a little practice to read the coefficients from the summary() output. For Example 2:
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 units (the intercept).
The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 365 + 37 = 402 units.
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 365 - 44 = 321 units.
The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 365 - 76 - 44 + 74 = 319 units.
As above, we can use the emmeans package to more easily see these coefficients:
emmeans(object = bestMod2, # your modelspecs =~ Cat2 + Cat3 + Cat2:Cat3, # your covariatestype ="response") # report coefficients on the response scale
So now we have a modelled value of our response for each level of our categorical covariate. For example:
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat1 is 365 +/- 25 units.
The modelled prediction for Resp2 when Cat2 is TypeB and Cat3 is Treat1 is 402 +/- 22 units.
The modelled prediction for Resp2 when Cat2 is TypeA and Cat3 is Treat2 is 321 +/- 27 units.
The modelled prediction for Resp2 when Cat2 is TypeC and Cat3 is Treat2 is 319 +/- 33 units.
Is each coefficient estimate different from zero?
Here we find out if there is an effect on Resp1 for all levels of Cat1 - are any of the coefficients statistically similar to 0?
forComp <-emmeans(object = bestMod2, # your modelspecs = pairwise ~ Cat2 + Cat3 + Cat2:Cat3, # your covariates, with the request for a pairwise testtype ="response") # report coefficients on the response scaletest(forComp)$emmeans # shows test of coefficients = 0
Here we see that the coefficients for all combinations of levels of our two covariates are significantly different than zero (t-test, P < 0.0001.)
Are the coefficient estimates different than one another across categorical levels?
As with Example 1, we can also find out which combinations of covariate levels are leading to statistically different model predictions in Resp2:
forComp <-emmeans(object = bestMod2, # your modelspecs = pairwise ~ Cat2 + Cat3 + Cat2:Cat3, # your covariates, with pairwise to indicate you want multiple comparisonstype ="response") # report coefficients on the response scaletest(forComp)
Again, the second table ($contrasts) shows the results of the multiple comparison testing. Based on a threshold p-value of 0.05, we can see that the effects at some combinations of our covariates are not statistically different from each other.
For example, the coefficient (predicted Resp2) when Cat3 = Treat1 and Cat2 = TypeA is 365 and this is 37 lower than the coefficient when Cat3 = Treat1 and Cat2 = TypeB, but this difference is not statistically different as p = 0.99 for this comparison.
On the other hand, some combinations of our covariates are statistically different from each other. For example, A comparison of modelled Resp2 when Cat2 = TypeD and Cat3 = Treat2 (245) vs. Cat2 = TypeB and Cat3 = Control (548) shows that they differ by 302 and are statistically different from one another with p < 0.001.
You can also get the results from emmeans visually with
plot(forComp,comparisons =TRUE)
Plotting modelled effects
You can visualize your results (model fit, uncertainty and observations) for Example 2 with:
# Set up your covariates for the visualized fitforCat2<-unique(myDat2$Cat2) # every level of your categorical covariateforCat3<-unique(myDat2$Cat3) # every level of your categorical covariate# create a data frame with your covariatesforVis<-expand.grid(Cat2 = forCat2, Cat3 = forCat3) # expand.grid() function makes sure you have all combinations of covariates# Get your model fit estimates at each value of your covariatesmodFit<-predict(bestMod2, # the modelnewdata = forVis, # the covariate valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data framelibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat2, # add observations to your plotmapping =aes(x = Cat2, y = Resp2, col = Cat3), position=position_jitterdodge(jitter.width=0.75, dodge.width=0.75)) +# control position of data points so they are easier to see on the plotgeom_errorbar(data = forVis, # add the uncertainty to your plotmapping =aes(x = Cat2, y = Fit, ymin = Lower, ymax = Upper, col = Cat3),position=position_dodge(width=0.75), # control position of data points so they are easier to see on the plotsize=1.2) +# control thickness of errorbar linegeom_point(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cat2, y = Fit, fill = Cat3), shape =22, # shape of pointsize=3, # size of pointcol ='black', # outline colour for pointposition=position_dodge(width=0.75)) +# control position of data points so they are easier to see on the plotylab("Resp2, (units)")+# change y-axis labelxlab("Cat2, (units)")+# change x-axis labellabs(fill="Cat3, (units)", col="Cat3, (units)") +# change legend titletheme_bw() # change ggplot theme
Reporting how well your model explains your response - Example 2
As for Example 1, you can report the explained deviance as the Likelihood Ratio \(R^2\):
So ~ 70% of the deviance in Resp2 is explained by Cat2 and Cat3 and the interaction between the two covariates.
Reporting covariate importance - Example 2
You can report the importance of each model term in explaining deviance in Resp2 with the sum of Akaike weights function (sw() in the MuMIn package), as we did in class:
Cat2 Cat3 Cat2:Cat3
Sum of weights: 1.00 1.00 0.96
N containing models: 3 3 1
This tells you that Cat2 and Cat3 appear in all models with any weight and that the interaction term Cat2:Cat3 is slightly less important in explaining deviance in Resp2, though it does appear in the best-specified model.
11.2.2 Predicting with Example 2
As with Example 1, it doesn’t make sense to make predictions with your model for Example 2 as you can’t make predictions to levels of your categorical covariates not already in your model.
11.3 Example 3: Resp3 ~ Cat4 + Cont5 + Cat4:Cont5
For Example #3, assume your best-specified model is that your response variable (Resp3) is explained by one categorical covariate (Cat4) and one continuous covariate (Cont5) as well as the interaction between the covariates:
Resp3 ~ Cat4 + Cont5 + Cat4:Cont5
Your best-specified model for example #3 is in an object called bestMod3:
For Example 3, you will report that your best-specified model is Resp3 ~ Cat4 + Cont5 + Cat4:Cont5. This model says that there is evidence that Cat4 and Cont5 explain variability in Resp3, and that there is an interaction between Cat4 and Cont5 - i.e. the effect of Cont5 on Resp3 depends on the value of Cat4. Remember that you can also use the output from dredgeOut3 to report evidence for how you picked the best-specified model, e.g. if you want to report more than one model.
The coefficients estimated for this model will include values of the response (Resp3) at each level of the categorical covariate (Cat4) as well as slopes describing the change in the response when the continuous covariate (Cont5) changes, and a different slope will be estimated for each level in Cat4 (this is because of the interaction in the model).
What are your coefficient estimates (with uncertainty)?
As above, you can take a look at your modelled coefficients with the summary() output:
The model prediction of Resp3 when Cat4 is Farm and Cont5 is 0 is 98.34 units56 (the intercept).
The model prediction of Resp3 when Cat4 is Urban and Cont5 is 0 is 98.34 - 0.24 = 98.1 units.
The model prediction of Resp3 when Cat4 is Wild and Cont5 is 0 is 98.34 - 0.87 = 97.5 units.
The slope of the relationship between Cont5 and Resp3 when Cat4 is Farm is -0.0029.57
The slope of the relationship between Cont5 and Resp3 when Cat4 is Urban is -0.0029 + 0.015 = 0.0121.
The slope of the relationship between Cont5 and Resp3 when Cat4 is Wild is -0.0029 + 0.031 = 0.0281.
Again, interpreting the coefficients from the summary() output is tedious and not necessary: You can use the emmeans package to give you the modelled response for each level of the categorical covariate (Cat4) directly:
emmeans(object = bestMod3, # your modelspecs =~ Cat4 + Cont5 + Cat4:Cont5, # your covariatestype ="response") # report coefficients on the response scale
Note that emmeans() sets our continuous covariate (Cont5) to the mean value of Cont5 (510 units). We can also control this with the at = argument:
emmeans(object = bestMod3, # your modelspecs =~ Cat4 + Cont5 + Cat4:Cont5, # your covariatestype ="response", # report coefficients on the response scaleat =list(Cont5 =0)) # control the value of your continuous covariate at which to make the coefficient estimates
By setting at = 0, you get the intercept - i.e. the modelled Resp3 when Cont5 = 0 for each level of Cat4, and this is what is reported in the summary() output.
Similarly, you can get the emmeans package to give you the slope coefficients for the continuous covariate (Cont5) using the emtrends() function, rather than interpretting them from the summary() output:
emtrends(bestMod, ~ Cat4, # your categorical covariatevar ="Cont5") # your continuous covariates
The emmeans() function will convert the coefficients (intercepts) from the link to the response scale. You can ask for this (as we did above) with the type = "response" argument. Note that it makes no difference for these examples as we are using a normal error distribution assumption and so the link and response scale are identical (i.e. when link = "identity").
In contrast, the emtrends() function does not convert the coefficients (slopes) to represent effects on the response scale. This is because emtrends() is reporting the slope of a straight line - the trend line on the link scale. But the line isn’t straight on the response scale.
We need to convert from the link to the response by hand when we use emtrends()^[remember, the coefficients for categorical covariates are also being converted from the link to the response scale. It is just that R is doing it for us in the emmeans() function when we add the argument `type = “response”). How we make the conversion, and interpret the result, will depend on our error distribution assumption:
Again, if we are using link = "identity" (as with a normal error distribution), the link and response scale are identical and no conversion is necessary. In this case, the coefficient tells you the absolute change in your response from a unit change in your covariate.
If you are using link = "log" (as with a poisson error distribution, and sometimes also a Gamma error distribution), you get your coefficient on the response scale by taking e to the coefficient on the link scale (with the R function exp()). This coefficient is called the rate ratio58 and it tells you the % change in the response for a unit change in your covariate.
But why do we convert coefficients from the link to response scale with ex when link = "log"?
If you are using link = "logit" (as with a binomial error distribution), you also get your coefficient on the response scale by taking e to the coefficient on the link scale (with the R function exp()), but here your coefficient on the response scale tells you the % change in the odds (probability of a success vs. the probability of a failure) given a unit change in your covariate - this estimate is called the odds ratio.
But why do we convert coefficients from the link to response scale with ex when link = "logit"?
Given a binomial model, \[
\begin{align}
log_e(\frac{p_i}{1-p_i}) &= \beta_1 \cdot Cov_i + \beta_0 \\
Resp_i &\sim binomial(p_i)
\end{align}
\]\(p_i\) is the probability of success and \(1-p_i\) is the probability of failure, and the odds are \(\frac{p_i}{1-p_i}\).
Then the odds ratio (OR), or % change in odds due to a change in covariate is:
For a model where link = "inverse", the interpretation on the response scale is less clear. In this case, you would use your model to make predictions and report the effects by describing these changes in predictions.
Is each coefficient estimate different from zero?
You can find out if the coefficients are different than zero, similar to above. You use emmeans() with test() when the covariate is categorical;
forInt <-emmeans(object = bestMod3, # your modelspecs =~ Cat4 + Cont5 + Cat4:Cont5, # your covariatestype ="response", # report coefficients on the response scaleat =list(Cont5 =0)) # control the value of your continuous covariate at which to make the coefficient estimatestest(forInt) # get test if coefficient is different than zero.
and we see that all coefficients are significantly different than zero (t-test, P < 0.0001).
For the continuous covariate, you use emtrends() with test():
forSlope <-emtrends(bestMod, # your model~ Cat4, # your categorical covariate, with a request for pairwise testingvar ="Cont5") # your continuous covariatetest(forSlope)
the slope associated with Cont5 when Cat4 = Farm is not different than 0 (p = 0.44). This means that there is no effect of Cont5 on Resp3 when Cat4 = Farm.
the slope associated with Cont5 when Cat4 = Urban is different than 0 (p = 0.0045) and positive (0.012). This means that Resp3 increases with Cont5 when Cat4 = Urban.
the slope associated with Cont5 when Cat4 = Wild is different than 0 (p < 0.0001) and positive (0.028). This means that Resp3 increases with Cont5 when Cat4 = Wild.
Are the coefficient estimates different than one another across categorical levels?
And as with the other examples, you can find which coefficients are significantly different from eachother across the levels of your categorical covariate. A rule here is to always check for differences in the coefficients associated with your continuous covariate(s) first as these differences are slope differences. If you have a difference in slope it is very likely you have intercept differences because the lines are not parallel.
You can also find out which slopes (i.e. the effect of Cont5 on Resp3) are different across the levels of Cat4 using emtrends() with a request for pairwise testing:
forCompSlope <-emtrends(bestMod, # your model pairwise ~ Cat4, # your categorical covariate, with a request for pairwise testingvar ="Cont5") # your continuous covariateforCompSlope
$emtrends
Cat4 Cont5.trend SE df lower.CL upper.CL
Farm -0.00293 0.00374 94 -0.01036 0.0045
Urban 0.01219 0.00418 94 0.00388 0.0205
Wild 0.02793 0.00492 94 0.01815 0.0377
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Farm - Urban -0.0151 0.00561 94 -2.694 0.0226
Farm - Wild -0.0309 0.00618 94 -4.991 <.0001
Urban - Wild -0.0157 0.00646 94 -2.437 0.0437
P value adjustment: tukey method for comparing a family of 3 estimates
From $contrasts you see that all slope estimates for all levels of Cat4 are significantly different from one another (0.0001 ≤ p ≤ 0.044).
Note that you can visualize these differences in slopes with:
plot(forCompSlope,comparisons =TRUE)
If slopes were similar, you will want to test if the levels of your categorical covariate (Cat4) have significantly different coefficients (intercepts) from each other with emmeans():
forCompInt <-emmeans(object = bestMod3, # your modelspecs = pairwise ~ Cat4 + Cont5 + Cat4:Cont5, # your covariates, with the request for a pairwise testtype ="response") # report coefficients on the response scaleforCompInt
From the bottom table ($contrasts), you can see that the intercept (modelled Resp2) modelled for each level in Cat4 is significantly different than every other level (all p-values are < 0.0001).
Plotting modelled effects
You can visualize your results (model fit, uncertainty and observations) for Example 3 with (the code is different than the other examples as you have a continuous covariate now):
# Set up your covariates for the visualized fitforCat4<-unique(myDat3$Cat4) # every level of your categorical covariateforCont5<-seq(from =min(myDat3$Cont5), to =max(myDat3$Cont5), by =1)# a sequence of numbers in your continuous covariate range# create a data frame with your covariatesforVis<-expand.grid(Cat4 = forCat4, Cont5 = forCont5) # expand.grid() function makes sure you have all combinations of covariates# Get your model fit estimates at each value of your covariatesmodFit<-predict(bestMod3, # the modelnewdata = forVis, # the covariate valuestype ="response", # make the predictions on the response scalese.fit =TRUE) # include uncertainty estimateforVis$Fit<-modFit$fit # add your fit to the data frameforVis$Upper<-modFit$fit+modFit$se.fit # add your uncertainty to the data frameforVis$Lower<-modFit$fit-modFit$se.fit # add your uncertainty to the data framelibrary(ggplot2) # load ggplot2 libraryggplot() +# start ggplotgeom_point(data = myDat3, # add observations to your plotmapping =aes(x = Cont5, y = Resp3, col = Cat4)) +# control position of data points so they are easier to see on the plotgeom_line(data = forVis, # add the modelled fit to your plotmapping =aes(x = Cont5, y = Fit, col = Cat4),size =1.2) +# control thickness of linegeom_line(data = forVis, # add uncertainty to your plot (upper line)mapping =aes(x = Cont5, y = Upper, col = Cat4),size =0.8, # control thickness of linelinetype =2) +# control style of linegeom_line(data = forVis, # add uncertainty to your plot (lower line)mapping =aes(x = Cont5, y = Lower, col = Cat4),size =0.8, # control thickness of linelinetype =2) +# control style of lineylab("Resp3, (units)") +# change y-axis labelxlab("Cont5, (units)") +# change x-axis labellabs(col="Cat4, (units)") +# change legend titletheme_bw() # change ggplot theme
As you saw in the previous section, the slope describing the effect of Cont5 on Resp3 when Cat4 is Farm was not significantly different than 0 (p = 0.44). Usually you would then not plot this line when plotting effects to avoid misinterpretation, but I’ll leave it on the plot for now.
Reporting how well your model explains your response - Example 3
As with the previous examples, you can report the explained deviance as the Likelihood Ratio \(R^2\):
r.squaredLR(bestMod3)
[1] 0.855657
attr(,"adj.r.squared")
[1] 0.8567834
So 86 % of the deviance in Resp3 is explained by variation in Cat4, Cont5 and the interaction between the two covariates.
Reporting covariate importance - Example 3
As with the previous examples, you can report the importance of each model term in explaining deviance in the response with the sum of Akaike weights function (sw() in the MuMIn package):
Cat4 Cont5 Cat4:Cont5
Sum of weights: 1 1 1
N containing models: 3 3 1
This tells you that Cat4, Cont5 appear in all models with any weight, indicating all terms have equal importance in explaining the deviance of the response.
11.3.2 Predicting with Example 3
Unlike with the other examples in this week’s notes, with Example 3, the presence of a continuous covariate (Cont5) means you can make predictions with your model - i.e. you can estimate what Resp3 might be if Cont5 was higher or lower. Because you also have a categorical covariate (Cat4) in your model, you do need to specify a level for Cat4 when you make your prediction.
Recall from Week 11 that we used the predict() function to make our predictions. For example, you can estimate the predicted Resp3 when Cont5 = 350 units and Cat4 = Urban with:
predict(object = bestMod, # your modelnewdata =data.frame(Cont5 =350, Cat4 ="Urban"), # the values of the covariates at which to make the predictiontype ="response", # to make the prediction on the response scale.se.fit =TRUE) # to include a measure of uncertainty around the prediction
So the predicted Resp3 when Cont5 = 350 units and Cat4 = Urban is 102.38 ± 0.82.
To get the predicted Resp3 when Cont5 = 1200 units when Cat4 is Wild, use:
predict(object = bestMod, # your modelnewdata =data.frame(Cont5 =1200, Cat4 ="Wild"), # the values of the covariates at which to make the predictiontype ="response", # to make the prediction on the response scale.se.fit =TRUE) # to include a measure of uncertainty around the prediction
you will be using the predict() function again in an upcoming chapter… when you predict!↩︎
note that the p-values reported in the coefficient table are the result of a test of whether the coefficient associated with Sp2 is different than that of Sp1 (p = 0.001), and if the effect of Sp3 is different than that of Sp1 (p = 0.017)), but a comparison of effects of Sp2 vs. Sp3 is missing. We will discuss this more further along in these notes.↩︎
note that the summary() output becomes easier to read if you force the starting model not to have an intercept with: startMod<-glm(formula = Resp1 ~ Cat1 - 1, data = myDat, family = gaussian(link="identity")). This is fine to do here where you only have categorical predictors but causes problems when you start having continuous predictors in your model as well. So we won’t be practicing this in class and instead will use the very flexible and helpful emmeans package to help us report coefficients from our statistical models.↩︎
“post hoc” is a Latin phrase meaning “after the event”↩︎
when the P-value is very low, R reports is as simple <.0001↩︎
note that now you get to assess the difference between coefficients for Sp2 and Sp3 which was missing in the summary() output above↩︎
The emmeans package is very flexible and has a lot of options as to how to make these corrections depending on your needs. Plenty more information is available here↩︎
also called the incidence rate ratio, incidence density ratio or relative risk↩︎
a reminder that the number e that the function exp() uses as its base is Euler’s number. It (along with the natural logarithm, log()) is used to model exponential growth or decay, and can be used here when you want models where the response can not be below zero.↩︎
remember, the coefficients for categorical predictors are also being converted from the link to the response scale. It is just that R is doing it for us in the emmeans() function when we add the argument `type = “response”)↩︎
i.e. effects associated with the continuous predictor↩︎
Notice I write “response(s)” in the title of this section - plural. It is possible to have multiple response variables and we’ll discuss this in class. For the focus of this course though, we’ll be working with one response variable↩︎
here, I’ll write “citation” where you would support your statements with the existing literature↩︎
As a reminder, in Week 9 we also discussed testing your hypothesis via p-values, as well as the limitations of this method when you have more than one covariate↩︎
AICc is the corrected Akaike Information Criterion which is more conservative than traditional AIC, i.e. models with more covariates need to increase the explained deviance quite a bit before the AICc metric improves↩︎
note that you can also include tables directly from R using the gt package↩︎
watch your number of significant units here. Make sure they make sense for the type of measurement. And be consistent↩︎
here equivalent to the traditional R2 as we have a normal error distribution assumption and linear shape assumption. The traditional R2 is found with 1-summary(bestMod)$deviance/summary(bestMod)$null.deviance↩︎
we’ll go over this in class today, but see also your course notes from week 9↩︎
note that the p-values reported in the coefficient table are the result of a test of whether the coefficient associated with Sp2 is different than that of Sp1 (p = 0.054), and if the effect of Sp3 is different than that of Sp1 (p = 0.0072), but a comparison of effects of Sp2 vs. Sp3 is missing. We’ll discuss this more further along in these notes.↩︎
note that the summary() output becomes easier to read if you force the starting model not to have an intercept with: startMod<-glm(formula = Resp1 ~ Cat1-1, data = myDat, family = gaussian(link="identity")). This is fine to do here where you only have categorical covariates but causes problems when you start having continuous covariates in your model as well. So we won’t be practicing this in class and instead will use the very flexible and helpful emmeans package to help us report coefficients from our statistical models.↩︎
here I say that your covariate “has an effect” on your response. Your covariate being in your best-specified model just means that there is evidence that variation in your covariate explains variation in your response. It’s up to you to think about the mechanisms underlying this relationship - i.e. is it causation or just correlation?↩︎
see section 11.6 in Crawley for more on multiple comparisons↩︎
when the p-value is very low, R reports is as simple <.0001↩︎
note that now we get to assess the difference between coefficients for Sp2 and Sp3 which was missing in the summary() output above↩︎
The emmeans package is very flexible and has a lot of options as to how to make these corrections depending on your needs. Plenty more information is available here↩︎